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Abstract 

This paper studies subordinate Ornstein-Uhlenbeck (OU) processes, i.e., OU diffusions 
time changed by Levy subordinators. We construct their sample path decomposition, show 
that they possess mean-reverting jumps, study their equivalent measure transformations, 
and the spectral representation of their transition semigroups in terms of Hcrmite expan- 
sions. As an application, we propose a new class of commodity models with mean-reverting 
jumps based on subordinate OU process. Further time changing by the integral of a CIR 
process plus a deterministic function of time, we induce stochastic volatility and time inho- 
mogcneity, such as seasonality, in the models. We obtain analytical solutions for commodity 
futures options in terms of Hcrmite expansions. The models are consistent with the ini- 
tial futures curve, exhibit Samuelson's maturity effect, and are flexible enough to capture a 
variety of implied volatility smile patterns observed in commodities futures options. 
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1 Introduction 



The contribution of this paper is two-fold. The first part studies subordinate Ornstein- 
Uhlenbeck (SubOU) processes. A SubOU process can be constructed by time changing an OU 
diffusion by a Levy subordinator. SubOU processes are Markov semimartingales with mean- 
reverting jumps. SubOU transition semigroups possess spectral representations in terms of 
Hermite expansions. As an application, the second part of the paper develops a new class 
of analytically tractable commodity models with mean-reverting jumps by modeling the com- 
modity spot price as the (scaled and compensated) exponential of a SubOU process. To model 
stochastic volatility and time inhomogeneity, such as seasonality, we further time change SubOU 
processes by the integral of the sum of an independent CIR diffusion and a deterministic func- 
tion of time. The resulting models have the following features: (1) mean-reverting jumps, (2) 
stochastic volatility, (3) time inhomogeneity, (4) analytical solutions for futures options in terms 
of Hermite expansions, (5) consistency with the initial futures curve, (6) Samuelson's maturity 
effect, and (6) flexibility to capture a variety of implied volatility smile patterns observed in 
commodity futures options. 

The mathematical part of the paper contains a self-contained presentation of SubOU pro- 
cesses. Section 2.1 defines SubOU semigroups as Bochner's subordinates of OU semigroups 
and gives explicit expressions for their infinitesimal generators based on the application of R.S. 
Phillips' theorem. This material is classical (see Schilling et al. (2010) for an excellent recent 
survey of Bochner's subordination and Albeverio and Rudiger (2003), (2005) for the treatment 
of SubOU semigroups in particular). Section 2.2 defines a class of SubOU Markov semimartin- 
gales, gives their local characteristics, proves uniqueness of the associated martingale problem, 
and proves the mean reversion property of their jumps. While the material in this section 
follows from the general semimartingale theory (our presentation follows Jacod and Shiryaev 
(2003)), it has not been presented in the literature in this form. Section 2.3 presents results 
on equivalent measure transformations for SubOU processes. In particular, a class of locally 
equivalent measure changes that transform one SubOU process into another SubOU process is 
characterized, along with a detailed treatment of some special cases important in applications. 
This section presents original results that, to the best of our knowledge, have not previously 
appeared in the literature. It serves as the basis for financial applications, characterizing equiv- 
alent martingale measures (EMMs) for this class of models. Section 2.4 presents the spectral 
decomposition of the SubOU semigroup in L 2 (M,m), where m is the Gaussian measure, in terms 
of Hermite expansions. The L 2 spectral theory of SubOU semigroups has been previously given 
by Albeverio and Rudiger (2003), (2005). We supplement it with pointwise convergence results 
and truncation error bounds for the expansion that are important for options pricing. 

The second part of the paper provides the development of our commodity futures model. 
Section 3.1 defines the model for the commodity spot price as the exponential of a SubOU 
process scaled and compensated so that, under Q, the mean spot price evolves along the fixed 
initial futures curve. We then explicitly solve for the futures dynamics under Q in the form of a 
martingale expansion with basis martingales associated with Hermite polynomials. Section 3.2 
demonstrates Samuelson's maturity effect in commodity futures in this class of models. Section 
3.3 derives explicit analytical solutions for futures options in terms of Hermite expansions. In 
section 4 we further time change SubOU processes to induce stochastic volatility and time 
inhomogeneity and study the resulting commodity futures models. In particular, we derive 
the futures price process, demonstrate Samuelson's maturity effect, and obtain solutions for 
futures options. In section 5 we discuss efficient model implementation based on recursions for 
Hermite polynomials and present model calibration examples to futures options on a variety 
of commodities, including metals, energies and agricultural . Appendix A contains a number 
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of results on the CIR process needed in the development of models with stochastic volatility. 
Proofs are collected in Appendix B. 

In the rest of this introduction we discuss relationships of models developed in this paper to 
the literature. We start with a brief survey of the commodity derivatives modeling literature. 
Mean reversion and jumps are two of the salient features of commodities prices (see monographs 



Eydeland and Wolyniec (2003), Geman (2005), and Geman (2008) for introduction to commod- 



ity and energy derivatives markets and modeling). Mean reversion in commodities markets 



is well documented in numerous empirical studies in the literature (e.g., Bessembinder et al 



(1995), Pindyck (2001), Casassus and Collin-Dufresne (2005)). To capture the mean reversion 



property, the classical commodity models are based on OU diffusions. The simplest such model 



is the exponential OU model of Schwartz (1997). In this model the commodity spot price is 



assumed to follow the exponential of an OU process with constant long-run mean level, rate of 
mean reversion, and volatility. While the OU process itself lives on the whole real line, taking 
the exponential leads to the positive process for the commodity spot price. The geometric OU 
model plays the same role in commodity markets that the geometric Brownian motion model 
plays in the equity markets, serving as the simplest analytically tractable commodity derivatives 
pricing model. Being inherently the spot price model, the futures curve is derived endogenously 
in this model and, hence, does not generally match the futures curve observed in the market. 



This situation is similar to the Vasicek (1977) model of the short interest rate, where the yield 



curve is derived endogenously in the model and does not generally match the market yield 
curve. Similar to how the Vasicek model is extended to match an arbitrary market yield curve 
by making the long-run mean level of the short rate time-dependent (e.g., Hull and White 
( |1993[ )), the exponential OU model can be extended to match an arbitrary market-observed 



futures curve (e.g., Clelow and Strickland (1999)). In this model futures prices of all maturities 
follow continuous martingales under Q. 

Along with mean reversion, discontinuous price movements (jumps) are another salient 
feature of commodity and energy markets. While jumps are a ubiquitous feature of all asset 
prices and financial variables, from equities to foreign exchange to interest rates, commodity 
and energy prices exhibit particularly large and frequent jumps, perhaps more so than other 
asset classes (see, e.g., Hilliard and Reis (1999), Deng ( |1999 ), Geman and Roncoroni (2005) 
for empirical evidence of jumps in commodity and energy prices). The question then arises as 
to how to extend commodity models based on mean-reverting OU diffusions to jumps. The 
first line of attack is to add a jump component to the diffusive mean-reverting component to 



form a jump-diffusion process similar to Merton (1976) classical jump-diffusion model widely 



used in equities. A variety of jump- diffusion models along these lines have been introduced in 



commodity markets (e.g., Hilliard and Reis (1998), Hilliard and Reis (1999), Deng (1999), Yan 



(2002), Benth and Saltyte Benth (2004), Geman and Roncoroni (2005), Andersen (2008) and 



Crosby (2008)). Virtually all of the jump-diffusion models in the literature, with the exception 



of Geman and Roncoroni (2005) and Andersen (2008), add state-independent jumps to the 



mean-reverting diffusion. The resulting models exhibit mean reversion due to the OU drift, 
but do not have mean reversion in their jump measure that remains state- independent. That 
is, upon arrival, the direction of the jump and the probability distribution of its amplitude are 
independent of the current state of the process. The drift acting upon the process between 
the jumps is forced to account for all of the mean reversion in these models. A model with 
mean reverting jumps would, in contrast, feature state-dependent mean reverting jumps with 
the jump direction and the jump amplitude dependent on the current state of the process. 



In contrast to jump- diffusion models with state-independent jumps, Geman and Roncoroni 



(2005) propose a jump-diffusion model with Poisson jumps independent of the diffusion com- 



ponent but with jump direction dependent on the pre-jump state of the process. They show 
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that such models capture some of the empirical properties of electricity price data. However, 
analytical solutions for futures options have not been obtained in their model. Andersen (2008) 
considers jump-diffusion processes with jumps driven by a continuous time Markov chain whose 
states are interpreted as different market regimes. Jumps in this framework are dependent on 
the regimes and, hence, are state dependent. However, option pricing in this regime-switching 
framework is generally highly non-trivial unless some simplifying assumptions are made. 

In this paper we take an alternative approach to the previous literature on commodity and 
energy models with jumps. Instead of adding state-independent jumps to the mean-reverting 
diffusion process, we time change the mean-reverting OU diffusion with a Levy subordinator to 
yield a pure jump or a jump-diffusion process (depending on whether or not the subordinator 
has a positive drift) with state-dependent and mean reverting jumps. As such, our models can 
be viewed as a commodity markets counterpart of the time-changed Levy process-based models 



in equity markets by Madan et al. (1998), Barndorff-Nielsen (1998), Geman et al. (2001), Carr 



et al. (2003), and Carr and Wu (2004). However, since mean reversion is the crucial feature 



of commodity markets, instead of time changing Brownian motion as in those references, we 
time change OU diffusions and obtain pure jump or jump-diffusion Markov semimartingales 
with state-dependent mean reverting jumps. Similar to Levy-based models in equity markets, 
our models based on SubOU processes calibrate well to a variety of implied volatility smiles in 
commodity markets when the maturity is fixed. 

To induce stochastic volatility (the need for stochastic volatility in energy markets has 
been advocated by Eydeland and Geman (1998])), we further time change these jump processes 
with the integral of an activity rate (stochastic volatility) that follows a CIR process. This 
is similar to the approach of Carr et al. (2003) and Carr and Wu (2004), but in contrast to 
those references we time change Markov jump processes that are generally not Levy processes. 
This yields pure jump or jump-diffusion models with stochastic volatility modulating jump 
amplitudes. To additionally introduce explicit time dependence to capture the term structure 
of at-the- money (ATM) volatilities observed in commodity futures options markets (e.g., the 
seasonality effects in volatility, as well as the sharply declining term structure of ATM volatility 
often seen in some commodity futures options), we add a purely deterministic function of time 
to the CIR activity rate (turning it into the so-called CIR+- 1- process, e.g., Brigo and Mercurio 



(2006)). Such models with mean-reverting jumps, stochastic volatility, and time dependence 
can be calibrated to the entire volatility surface across both the strike and maturity dimensions. 
To conclude this introduction, we discuss analytical and computational aspects. While the 



time changed Levy models of Carr et al. 



of Fourier analysis (see Carr and Madan 



(2003) lead to fast and efficient option pricing by means 



1999) for the Fast Fourier Transform methodology and 



Feng and Linetsky| ( 2008 ) and Feng and Linetsky ( 2009 ) for the closely related Hilbert trans- 



form methodology), our time changed OU models also lead to analytical option pricing, but by 
different mathematical means. While in the context of Levy processes one exploits the explicit 
knowledge of the characteristic function, in the context of OU processes we exploit the explicit 
knowledge of the eigenfunction expansion of the SubOU transition semigroup. The eigenfunc- 
tion expansion method is a powerful tool for pricing contingent claims written on symmetric 
Markov processes (see Linetsky (2004) and Linetsky (2007) for surveys). It is particularly well 
suited to time changes since the time variable enters the eigenfunction expansion of the tran- 
sition semigroup only through the exponentials e~ Xnt and, after the time change with a Levy 
subordinator, the eigenfunction expansion has the same form as for the original process, but 



with e Xnt replaced with e 



-0(A„)t 



where is the Laplace exponent of the subordinator. We 



note that the seminal paper by Bochner ( |1949 ) already contained this observation (see Eq.(ll) 



m 



Bochner (1949); further see Albeverio and Rudiger (2003), (2005) for the mathematical de- 



velopment of subordination of symmetric Markov processes). In Mathematical Finance, this 
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observation has been previously exploited by Albane se and Kuznetsov ( 2004|) in the context of 
volatility smile modeling for equities, by Boyarchenko and Levendorskii (2007) in the context 
of interest rate modeling, and by Mendoza et al. (2010) in the context of unified credit-equity 
modeling. 



2 Subordinate Ornstein-Uhlenbeck Processes 



2.1 SubOU Semigroups 

We start with an OU semigroup (Vt)t>o defined on !Bb(IR) (the space of bounded Borel 
measurable functions), where Vtf{x) = J K f(y)p(t,x,y)dy with the OU transition kernel: 



K \ 



cxp 



(y - x + (z - 0)(1 - e~ Kt ))' 



-2Kt\ 



£l(l 



-2Kt\ 



}■ 



(2.1) 



p(t, x, y) is the transition density of an OU diffusion with the rate of mean reversion n > 0, long- 
run level 9 6 M, and volatility a > 0. (Vt)t>o is a strongly continuous contraction semigroup 
on QS&QR). Restricted to Cq(M) (the space of continuous functions vanishing at infinity), it 
is a Feller semigroup, and C£°(R) is a core of the domain D(Q) of its infinitesimal generator 
Q acting on C|(R) (the subsc ript c stands for fu nctions with compact support) by Qfix) = 
k{6 - x)f(x) + \a 2 f"{x) (c.f. |Duffie et~al| ( |2003[ ) Theorem 2.7). 

Consider a vaguely continuous convolution semigroup (qt)t>o of probability measures on 
R+ (c.f. Schilling et al. (2010) Definition 5.1). For each t, (ft([0,oo)) = 1 (we consider only 



conservative case in this paper), and its Laplace transform is given by the Levy-Khint chine 
formula with the Laplace exponent (j>(\) defined for all A > 0: 



-As 



[0,oo) 



q t (ds) 



-t<f>{A) 



0(A) = 7 A + 



[0,oo) 



(1 



)v(ds) 



with drift 7^0 and Levy measure v satisfying the integrability condition Jj oo ^(sAl)u(ds) < oo. 
(qt)t>o is the family of transition probabilities of a subordinator, i.e., a non-negative Levy process 



starting at the origin (c.f. Bertoin (1996) or Schilling et al. (2010)). 

We define a subordinate semigroup {Vf)t>o on ?Bfe(]R) as the Bochner integral: 



Vff(x) 



f V s f{x)q t {ds) 

J[0,oo) 



This procedure is cal 
From Schilling et al. 



e d Boc hner's subordination (c.f. Schilling et al. ( 2010[ ) Definition 12.2). 
(2010) Proposition 12.1, the subordinate semigroup ivf)t>0 is also a 
strongly continuous contraction semigroup on Q3b(R). We call it the SubOU semigroup with 
generating tuple (K,9,a,"f,u). The superscript </> in (Vf)t>o signifies that it is constructed by 
subordinating the semigroup (Vt)t>o with the convolution semigroup of a subordinator with the 
Laplace exponent 4>. 



From Jacob (2001) Corollary 4.3.4, a Feller semigroup remains a Feller semigroup after 
subordination. It implies that (Vf)t>o restricted to Co(M) is Feller. Its infinitesimal generator 
is given by Phillips' Theorem ( |Sato (1999) Theorem 32.1). The assertion on its core comes from 
Sato (1999) Proposition 32.5 (ii) and the fact that C£°(M) is a core of D{Q). We summarize 



these results in the following. 
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Theorem 2.1. (i) A SubOU semigroup with generating tuple (k, 9, a, 7, v) is a Feller semigroup. 
(ii)Let be its infinitesimal generator. Then C™(R) is a core of D(Q^), C 2 {R) C and 
for any f G C C 2 (M), 

Q*f&) = \l° 2 f'\x) + Hx)f'(x) + / (f(x + y)- f(x) - yl M<1} f'(x)) U(x, dy), 
with the state- dependent Levy measure II(x, dy) = 7r(x, y)dy with density defined for all y 7^ 



(2.2) 



n(x,y)=l p(s,x,x + y)v{ds), 
'[0,00) 



where p(t,x,y) is the OU transition density (2.1). The drift with respect to the truncation 
function yl-n^i^u is 



b(x) = 7k(6> -x) + 



[0,oo) V-/{[v|<l} 



yp(s,x,x + y)dy v(ds). 



Remark 2.1. (i) On C?(R) can be represented as a pseudo-differential operator (PDO) (see 
for example |Schnurr| |2009[) Corollary 1.21) G*f{x) = -p(x,D)f(x) = - hpfaOf 



where /(£) = ^- L e~ l ^ x f(x)dx is the Fourier transform of f(x), and p(x,£) is called the symbol 
of the PDO and is expressed as p(x, £) = \^o 2 ^ — ib(x)^ — j y , Q {e^ v — 1 — n(x, dy). 

Note that p(x,£) is a continuous negative definite function (CNDF) for each x (c.f. Jacob ([2001) 
Definition 3.6.5). 

(ii) n(x, y) satisfies the condition J y _^ (y 2 A l)ir(x,y)dy < 00 for each x. This is a direct result 
from the representation theorem for CNDF. See Jacob (2001) Theorem 3.7.7. 



(iii) The Levy measure of the SubOU semigroup has finite activity if and only if the Levy 
measure of the subordinator has finite activity, which is justified by interchanging the order of 
integration in f y ^ Q J, Q s p(s, x, x + y)v(ds)dy by Tonelli's Theorem. 

(iv) In general, it is not true that we can interchange the order of integration in 
J[o 00) f\y\<i yp( s > x ,x + y)dyv{ds). However, if the Levy density satisfies the integrability condi- 
tion fiyi<i \y\ 7T (x,y)dy < 00, then the interchange is valid and the truncation is not needed. It 
can be shown that this integrability condition is equivalent to y/su(ds) < 00 (if x 7^ 9) and 
Jq v(ds) < 00 (if x = 9). In this case the generator takes the simpler form on C 2 (M): 

G*f(x) = l -ia 2 f"(x) + 1 k{9 - x)f'{x) + [ (f(x + y)- f{x)) ir(x, y)dy. 



2.2 SubOU Processes as Markov Semimaringales 

Definition 2.1. A time-homogeneous Markov process (f2, J 7 , (J-t)t>o, X, ¥ x ) x ^ with state space 
(R, is called a subordinate OU (SubOU) process with generating tuple 

(k,9,ct,7,v) if its semigroup is a SubOU semigroup with the same generating tuple. 

Since a SubOU semigroup is Feller, a SubOU process is a Feller process. Every Feller process 



has a cadlag modification (c.f. Jacob (2005) Theorem 3.4.9 or Revuz and Yor (1999) Theorem 



III. 2. 7), so immediately we have the following 

Corollary 2.1. Every SubOU process admits a cadlag modification. 
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We will always consider cadlag SubOU processes in this paper. From now on, without 
explicit mention, we will assume that (A, (P x ) x eM.) is the canonical realization of a given SubOU 
semigroup defined on ( »,J 70 , (-F t °)t>o), where » = D(R) (the Skorohod space of cadlag functions 
with values in R, c.f. Jacod and Shiryaev (2003) Definition VI. 1.1), J7 = a(X s ,s < t), and 



| Schnurr (2009) gives an excellent discussion on the connection between cadlag Feller pro- 
cesses and semimartingales. The Feller property of the SubOU process together with C^°(R) C 
D(Q^) allows us to conclude that it is a semimartingale w.r.t. every P x with i£M (c.f. Schnurr 



(2009) Theorem 3.1). From Schnurr (2009) Theorem 3.14, the pseudo-differential operator rep- 
resentation of the infinitesimal generator gives us the triplet (B, C, II) of semimartingale 
characteristics of the SubOU process. For the definition of semimartingale characteristics see 



Jacod and Shiryaev ( 2003 ) Chapter II 



Theorem 2.2. (i) The SubOU process (O, J 70 , {J^)t>o, X, P^eR with generating tuple 
(k, 6, cr, 7, v) is a semimartingale w.r.t. every P x and admits semimartingale characteristics 



(B,C, II) w.r.t to the truncation function h(x) 



xl {N<i}' where 



B t (u) 



o 



jk(9 - X s .(oj)) + 



o J{\y\<i} 



yp(u; A s _ (to), X s - (cj) + y)dyv{du) 



ds, 



C t (ui) = ja 2 t, U(u,dt,dy) = ir(X t -(uj),y)dtdy, 



where n(x,y) is given in Theorem 2.1. 

(ii) Denote by ji x the integer-valued random measure associated with the jumps of X (c.f. 



Jacod and Shiryaev (200S\ ) Proposition II. 1.1 6) and X c the continuous local martingale part of 
X. Then X has the following sample path decomposition (under the starting point x): 



X t {oo) = x + B t (uj)+ X£(cj) + h(x) * (/i x - n) t (w) + (x - h{x)) * (w) 



(2.3) 



with the quadratic variation of the continuous part [X c , X c }t(ui) = Ct(to) = ja 2 t (* denotes 
integration w.r.t. a random measure). 

(Hi) If X' is an M-valued semimartingale defined on some filtered probability space 
(£l',J-', (Jj)t>o,P') withP' (X' = x) = 1 and with the semimartingale characteristics (B',C, IT) 
given in (1), where X is replaced by X' , then P' o X'^ 1 = P x . 



The proof of Part (iii) of Theorem 2.2 is given in Appendix B. It essentially says the solution 



to the Martingale Problem in the canonical space setting as defined in Jacod and Shiryaev ( 2003 ) 



Definition III. 2.4 is unique. This is a key result for the study of locally equivalent measure 



changes for SubOU processes in section 2.3 



Remark 2.2. (i) It is clear that the SubOU process is a jump-diffusion process if 7 > and a 
pure jump process if 7 = 0. 

(ii) If Jq 1 yfsv{ds) < 00, then Jj 2/ i <1 \y\'x(x,y)dy < 00 for all x 7^ 9. Hence \h(x)\ * II G ^u, c i 
which implies h{x)*{^ x — II) = h{x) *[i x — h(x) *II (c.f. Jacod and Shiryaev (2003) Proposition 
II.1.28) and 

X t {oj) = x+ / jK(0-X s -.(u))ds + X£(u) + x*[m x (cu). 



Hence, in this case, the jump part of the SubOU process is of finite variation. 

A SubOU process is a process with mean- reverting jumps. The mean reversion property of 
the state-dependent SubOU Levy measure Tl(x, •) is characterized in the following. 
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Theorem 2.3. For any y > 0, we have 

(i) If x > 9, then tt(x, —y) > ir(x, y), and Ii{x, (— oo, —y)) > IL{x, (y, oo)). 

(ii) If x < 9, then tt(x, —y) < ir(x, y), and Ii{x, (— oo, —y)) < H{x, (y, oo)). 
(in) If x = 9, then tt(x, —y) = ir(x, y), and Ii{x, (— oo, —y)) = n(x, (y, oo)). 

This theorem tells us that when the current state x is above (below) the long-run level 9, a 
downward (upward) jump is more likely to occur. When x = 9, the intensity of downward and 
upward jumps are equal. This mean-reverting nature of jumps makes SubOU processes a natural 
candidate for modeling mean-reverting prices and other financial variables. If 7 = 0, a SubOU 
process is a pure jump process with mean-reverting jumps. If 7 > 0, it is a jump- diffusion 
process with mean-reverting diffusion drift and mean-reverting jumps. 

Figure [T] plots SubOU Levy densities when v are Levy measures of a compound Poisson 
process with exponential jump sizes and an inverse Gaussian (IG) process. 

Remark 2.3. Time Change Interpretation of Bochner's Subordination The semigroup 
{Vt)t>0 gives rise to an OU diffusion process X. The vaguely continuous convolution semigroup 
of probability measures (qt)t>o gives rise to a subordinator T. Assume that both X and T 
are defined on the same probability space and are independent. Then the time changed or 
subordinate process xf := Xx t is again a Markov process. By independence of X and T, the 
associated operator semigroup is given by 

V?f(x)=E[f(X Tt )}= [ E x [f(X 8 )]q t (ds) = [ V s f(x)q t (ds). 

J [0,0a) J [0,00) 

That is, xf is a SubOU process according to our definition, and Bochner's subordination can 
be interpreted as a stochastic time change with respect to an independent subordinator (cf. 
Schilling, Song and Vondracek (2010) p. 141). 

Remark 2.4. SubOU Markov semimartingales admit a representation in terms of a Brownian 
motion and an independent Poisson random measure. Explicit expressions follow from Cinlar 
and Jacod (1981) Theorem 3.13 and are omitted due to space constraints. 



2.3 Equivalent Measure Transformations for SubOU Processes 



For building financial models based on SubOU processes, we are interested in locally equiv- 
alent measure changes^ that transform a SubOU process with a given generating tuple into 
another SubOU process with another generating tuple. We can then build financial models 
with SubOU processes under both the physical and the risk- neutral measures, and determine 
how the generating tuple of the SubOU process changes under the measure change. 

As before, 0, is the space of all cadlag functions taking values in R. In this section we follow 



Jacod and Shiryaev ( |2003[ ). In order to use their results, we use the right-continuous version of 

J 70 . Let X be the canonical process. 



the filtration (Ft)t>o with Ft 



Ft + and F 



Vt>o ?t 



It is clear that if X is a SubOU process, it is also Markov and a SubOU process w.r.t. (F)t>o- 
We fix the truncation function h(x) = xl| x< 1 j. Let Pq be a probability measure on (fi, Fq) 



taken to be the initial distribution. Following Jacod and Shiryaev (2003) Definition III. 2. 4, we 
call a probability measure P on (fi, F, (Ft)t>o) a solution to the martingale problem associated 
with (Fo,X) and (¥q;B,C, u), where (B,C,v) are given semimartingale characteristics, if the 

Two probability measures Pi and P2 on a filtered probability space (£1, (J-t)t>o) are said to be locally 
equivalent, if Pi|j^ ~ P2l.Fi for each t > 0, where P|jr t is the restriction of measure P on the cr-filed T t - 
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Figure 1: State-dependent Levy densities of SubOU processes with 9 = 0.2, k = 1, and a = 0.6 
when v is the Levy measure of a compound Poisson process with exponential jumps (arrival 
rate a = 2, reciprocal of mean jump size rj = 1) and an inverse Gaussian process (mean rate 
fi = 1, variance rate v = 1) with x = —1, 0.2, 1. To emphasize the value of the current state, 
the horizontal axis plots the post jump state (not the jump size). 
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following hold: (i) the restriction P|^ o = Po; (p) X is a semimartingale on the stochastic 

basis (0, J 7 , (J'jjoo, IP) with characteristics (B,C,u) relative to the truncation function h. The 
following proposition is crucial in proving the necessary and sufficient conditions for locally 
equivalent measure change. 

Proposition 2.1. Let (B,C,H) be the SubOU semimartingale characteristics defined in The- 
2.2. The solution to the martingale problem (<r(Xo), X|Pq; -B, C, II) exists and is unique. 



orem 



Moreover, local uniqueness holds. 



See Jacod and Shiryaev (2003) Definition III. 2. 35 for the definition of local uniqueness. The 
existence of the solution is quite obvious. Given a SubOU semigroup with generating tuple 
corresponding to the given SubOU semimartingale characteristics (B,C,v), we can construct 
a time-homogeneous universal Markov process on the space of cadlag functions taking values 
in R. Such a process is a semimartingale with characteristics (B, C, v) by Theorem 2.2 under 
every P x , and set ¥(A) = j¥ x (A)¥o(dx) for any A G T . The uniqueness follows from part (iii) 
of Theorem |2.2| The local uniqueness is a result of uniqueness and the Markov property of the 
process by Jacod and Shiryaev (2003) Theorem III. 2. 40. We then have the following. 



Theorem 2.4. Let P and P' be two probability measures on (£1, J-, (J-t)t>o) such that the canon- 
ical process is a SubOU process with generating tuples (K,9,a,j,v) and (k',0 1 , a' ,7' ,u'), respec- 
tively, and with initial distributions Po and P 0; respectively. Then the following two statements 
are equivalent. 

(1) P and P' are locally equivalent, i.e., ¥\p t ~ ^'\p t f or ever U t >0. 



(2) The following conditions are satisfied: 

(1) P ~ P' ; (ii) 7d 2 = 7V 2 ; 2 

(iii) For every the Hellinger condition f , Q ^tt'(x, y) — ^ir(x, y)) dy < cc holds, 

where tt(x, •) and ir'(x,-) are defined as in Theorem\2.l\ 



(3) Furthermore, suppose these conditions are satisfied. Define 
{in'9 1 - 7K#) - (7V - yK)X t -(u) 



AH 



1 {7 ^0}, and Y(u,t,y) r- 



70- ir{X t -(u}),y)' 

Let X c and \x x denote the continuous martingale part of X and the jump measure associated 
with X . Then N = (3 ■ X c + (Y — 1) * (fi X — II) is a F '-local martingale, and the Radon-Nikodym 
density process D o/P' w.r.t. P equals to the Doleans-Dade stochastic exponential S{N) of N. 

Remark 2.5. Define ip(x,y) := In (V(x, y)/ir(x, y)). The Hellinger condition L^q {y/^'ix, y) — 



Sato 



(1999)): 



{y:tp(x,y)<-l} 



ir(x, y)dy < 00. 



i/7r(x, y)) 2 dy < 00 is equivalent to the following (similar to Remark 33.3 of 
/ ¥(x,y) 2 n(x,y)dy < 00, / n'(x, y)dy < 00, 

Intuitively, when tt and tt' both have infinite activity, the Hellinger condition says that the 
region where large perturbations of the jump density occurs should not be arbitrarily close to 
the origin. 

Remark 2.6. The limiting case k = corresponds to the subordinate Brownian motion without 
drift. Theorem 2.4 still holds when k = and/or k' = 0. When k > and k' = 0, it 
characterizes locally equivalent measure transformations of SubOU processes into subordinate 
Brownian motions without drift. When k = k' = 0, Theorem 2.4 reduce to the special case of 
Theorem 33.1 of Sato ( 1999 ) for Levy processes specialized to the case of subordinate Brownian 
motions. 
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The general Hellinger condition is difficult to check. We wish to derive restrictions it places 
on SubOU generating tuples that can be transformed into each other under locally equivalent 
measure changes. It can be easily shown that: (i) If both v and v' are Levy measures of finite 
activity subordinators, then the Hellinger condition is automatically satisfied, (ii) If v is a 
Levy measure of a finite activity subordinator and v' is a Levy measure of an infinite activity 
subordinator (or vice versa), then the Hellinger condition is not satisfied. Thus, equivalent 
measure changes cannot transform a SubOU process with a finite activity subordinator into a 
SubOU process with an infinite activity subordinator, and vice versa. 

We now investigate the case when v and v' are Levy measures of infinite activity subordina- 
tors. To verify the Hellinger condition in this case, we need to study the asymptotic behavior 
of the SubOU Levy density Tr(x,y) given in Eq.(2.2) as y — > 0. The following proposition 
shows that it is equivalent to the asymptotic behavior of the Levy density of some subordinated 
Brownian motion. 

Proposition 2.2. Let Tr(x,y) be the Levy density of a SubOU process with generating tuple 
(k, 6, a, 7, is). Suppose ir(x,y) — > oo as y — > 0. For each fixed x £ M., let Tf x (y) be the Levy 
density of a subordinate Brownian motion starting at with drift k{9 — x), volatility a, and the 
same 7 and v. Then lim^o tt(x, y)/n x {y) = 1- 

We can further show that the asymptotics of the Levy density of subordinate Brownian 
motion does not depend on drift. We then have. 

Proposition 2.3. The asymptotics ofir(x,y) as y — > does not depend on k, 9, and x. 

Hence, k and 9 can be freely changed by locally equivalent measure changes. In particular, 
k > can be changed to k = by a locally equivalent measure change. The problem of 
investigating the Hellinger condition now reduces to finding the asymptotics of the Levy density 



of a subordinate Brownian motion. Song and Vondracek (20091 is an excellent reference on the 



potential theory of subordinate Brownian motions and provides many examples of subordinators 
and asymptotics of the Levy densities of subordinate Brownian motions. If the Levy measure 
v of the subordinator has a density v(s), then, in general, we have Proposition 2.4 to compute 
the asymptotics of the Levy density of the subordinate Brownian motion 

1 



7r(y) :-- 



y 

2a 2 



iro~ 2 s 



v{s)ds 



(2.4) 



asy 



0. Proposition 2.4 gives the asymptotics under two different types of sufficient conditions. 
The first sufficient condition is based on Lemma 3.3 of Song and Vondracek (2009). The 



applicability of their Lemma 3.3 is not restricted to Levy densities of subordinators. However, 
in this case some of their conditions are not necessary. Below we give a more general result for 
this case. The second sufficient condition is a restriction of the Levy density of the subordinator 
to the class of completely monotone function^] (see, for example Schilling et al. (2010) for its 
characterization and properties). This result is proved in Theorem 2.6 in Kim et al. (2010). 



Proposition 2.4. Let z/(s) be the Levy density of a subordinator. Suppose there exist constants 



Co > and \ < (3 < 2 and a function I : (0, 00) 



MS) 



(0, 00) slowly varying at infinity such that 
0. (2.5) 



as s 



Let n(y) be defined as in (2.4). Suppose one of the following two conditions is satisfied: 

2 A completely monotone function / : (0, 00) 1— > R is a C°° function such that (~l)"f^(x) > for n 
0,1,2,.-.. 



A function I defined in a neighborhood of infinity is called slowly varying at infinity if lim^ 



all a > 0. 



t(ax) 



1 for 
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(1) Let g : (0, oo) —> (0, oo) be a function such that f. 



oo J. 



3 



s g(s)ds < oo. Assume there is 



also some £ > such that fg^(y,s) < g(s) for all y, s > 0, where the auxiliary function 

£(-) 

fl^{y,s) is defined by fi^{y,s) := 2 J 2s if y < § and otherwise for any function t 



slowly varying at infinity and any £ > 0. 
(%) i>(s) is a completely monotone function. 
Then 



coT(P 



7r(2cr 2 ) 1 -^ | y p/5-l£(^) 



as y — > 0. 



Remark 2.7. For slowly varying functions and regularly varying functions see Bingham et al. 

l 



(1987). Every regularly varying functior 
number /3 and 



at zero can be written in the form 



TO) 



for some real 



slowly varying at infinity (c.f. Bingham et al. (1987) Theorem 1.4.1). Hence 



the assumption on the asymptotics (2.5) is very general. Also note that from Bingham et al. 

l 



1987) Proposition 1.3.6, 



TO) 



oo as s — t- 0, so we are dealing with subordinators whose 



evy density tends to infinity at 0. 



For a subordinator with Levy density, if (2.5) is satisfied, there is a close connection be- 



tween the Blumenthal-Getoor (BG) index and the parameter (3 in Proposition 2.4 when fj > 1. 
For any subordinator with Levy measure v{ds) its BG index is defined by p := inf{a > : 
J| s |<i s a "(ds) < oo}. 

Proposition 2.5. (1) Suppose (2.5) holds with f3 > 1. Then the BG index is equal to f3 



1. 



(2) Suppose the conditions in Proposition 2.4 are satisfied for two subordinators with /3 > 1 and 
j3' > 1. Then the Hettinger condition implies their BG indexes are equal. 



We now apply Proposition 2.4 to the key example important in financial applications. 



Example 2.1. Tempered Stable Subordinators. Consider the tempered stable family of 
Levy measures u{s) = Cs~ 1 ~ p e~ vs , where C > 0, p < 1, n > 0. The limiting stable family has 
rj = and p G (0, 1). The tempered stable cases with p > (p < 0) give rise to subordinators 



p = ( 


Madan et al. 


Nielsen\ (1998 


)), and 



(Barndorff- 



Nielsen\ (1998)), and the compound Poisson subordinator with exponential jumps with p 
and rj > 0. For this family, the Laplace exponent is given by the following. 

f jX-Cr{-p)[(X + rj)P-rfP], p^0 
7 A + Cln(l + A/77), p = 0' 



0(A) 



-1 



(2.6) 



where T(-) is the Gamma function. 

For tempered stable subordinators with p > 
holds with co = C, j3 = 1 + p, t{x) = 1, g(s) = 



2., 



condition (1) 



-g; it is clear that Proposition 
1, and!; chosen arbitrarily. Condition (2) also 
holds because the Levy density of the subordinator is completely monotone. Hence we have 



7r(y) 



CT{p+\)(2a 2 Y 



ir\y 



2p+l 



as y 



0. 



From Proposition 2.2, ir(x,y) has the same asymptotics. It is now straightforward to show that 
Theorem 2.4 reduces to the following result for Sub OU processes with tempered stable subordi- 
nators with drift. 



4 A function / defined in a neighborhood of infinity is called regularly varying at infinity with index p G 
lim^-^oo ^f(*) = ^ P f° r au ^ > 0. It is called regularly varying at if /(j) is regularly varying at 00. 
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Corollary 2.2. Consider the setting in Theorem \2.J\ Suppose v and v' belong to the tempered 
stable family with parameters (C,p,rj) and (C',p',r]') with p,p' > 0, respectively. Then P|_^ ~ 
P / 1 ~ for every t > if and only if Pq ~ Pq an< ^ ^ e following equalities hold: 



.2 _ „,/_/2 p = ^ Ca 2p = c / a /2 P _ 



7cr =70" 



Thus, if we have SubOU processes with tempered stable subordinators with drift under 
both the physical and the risk-neutral measure, the p parameter p must remain the same under 
both measures, C and C are related by C = C(o~/a') 2p , the subordinator drifts and the OU 
volatilities are related by 7V 2 = 70" 2 , and the OU drift parameters 9 and k and 9' and k! can 
be arbitrarily changed. 

For other examples of Levy densities, where, e.g., £(x) = (ln(l + x)) a , one can also use 



Proposition 2.4 Sec Song and Vondracek (2009) section 2 for examples of subordinators and 



section 3 for the asymptotics of the Levy density of subordinate Brownian motions. Replace 
4 in their formulas by 2cr 2 to coincide with our notation. Once the asymptotics of the Levy 
density is determined, the Hellinger condition can be reduced to a simple relationship for the 
parameters similar to Corollary 2.2 for SubOU processes with tempered stable Levy densities. 

We are also interested in the following question: if under some measure P the semimartin- 
gale X is a SubOU process with generating tuple (k, 9, a, 7, u), characterize all measures P' 
locally equivalent to P. In particular, we are interested in conditions on the semimartingale 
characteristics of X under P'. The following result answers this question. 

Theorem 2.5. Let P and P' be two probability measures on (fi, J 7 , (J-~t)t>o) with initial distri- 
butions Po and Pg, respectively. Suppose under P, the canonical process X is a SubOU process 
with generating tuple (k, 9, a, 7, v) and local characteristics (B,C, II). Suppose under P', X is 
a semimartingale with local characteristics (B' ,C' ,H'). IfF' andF are locally equivalent, then 
there exists a nonnegative predictable function Y(lj, t, x) and a predictable process /3 such that: 

B' t (u) = B t (u)) +7(j 2 / f3 s (u)ds+ / yl{\ y \< ;1 y(Y(s,u,y) - l)Tr(X s _(u),y)dyds, 

JO J[0,t]xR 

C' t {uj) = -ya 2 t, W(u,ds,dy) = Y(u,s,y)U(uj,ds,dy), (2.7) 
,-t ,-t 

Jo Jo 

t 

Jy+0 



I \(3 s (u))\ds < 00, / (3;(cu)ds < 00, (2.8) 
Jo 

|y 1 {|y|<i}( y ( s '^'y) ~ 1 )\ 7T ( x s-{u),y)dyds < 00, (2.9) 



(a/Y(s, lo, y) — l) 2 n(X s -(u),y)dyds < 00, (Hellinger condition) (2-10) 

Jy^o 

P' and F-a.s. for all t > 0. Define N = f3 • X c + (Y - 1) * (fi x - II). Then the density process 
Z o/P' w.r.t. P is the Doleans-Dade stochastic exponential <§{N) of N. 



2.4 The Spectral Representation of the SubOU Semigroup 

The OU and SubOU processes are stationary with the Gaussian stationary density 

/ K «(A-*) 2 
m X = \/ 2 e ^ • 

Consider the Hilbert space L 2 (M,m) with the inner product (/,<?) = L f(x)g(x)m(x)dx, and 
denote by || • || the L 2 -norm. The OU and SubOU semigroups are both symmetric semigroups 
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in L 2 (M,m), i.e. (V t f,g) = (f,V t g) and (V?f,g) = (f,Vfg) for any f,g G L 2 (R,m). Their 
spectral decompositions in L 2 (R,m) are available in closed form. 

Theorem 2.6. (1) The OU semigroup has the following eigenfunction expansion in L 2 (M,m): 



V t f(x) = J2e~ Knt fnMx), f^L 2 



,m 



t > 0, 



n=0 



1 



(x-9) 



n = 0,l, 



\/2 n n! \ cr 
and expansion coefficients f n = (f,ip n ). 

(2) The SubOU semigroup has the following eigenfunction expansion in L 2 (IR,m): 



(2.11) 



with the orthonormal eigenf unctions expressed in terms of Hermite polynomials (see, e.g., Lebe- 
dev (1965)) 



(2.12) 



V?f(x) = f2 e ~ HKn)t fnMx), /G£ 2 (M,m), t >0 



(2.13) 



n=0 



with the same eigenf unctions and expansion coefficients as the OU semigroup. 



General results for the spectral representation of one-dimensional diffusions go back to the 
fundamental work of McKean (1956). For each t, the OU transition semigroup operator Vt 
has a purely discrete spectrum with eigenvalues {e~ Knt ,n = 0, 1, ...}. The explicit form of the 
eigenfunction expansion of the OU semigroup in terms of Hermite polynomials is well known and 



can be found in many references 


, including Wong 


(1964 


, Karlin and Taylor (1981), Schoutens 


(2000 


), Bakry and Mazet 


(2004 


), 


Alberverio and Riidiger 


(2003 


), 


Alberverio and Riidiger 


(2005 


), 



and Gorovoi and Linetsky (2004) p. 62. The general spectral representation of the transition 



semigroup of a symmetric Markov process can be found in Fukushima et al. (1994). Bochner 
subordination replaces the eigenvalues e~ An< with e~^ Xn ^, where (j> is the Laplace exponent 
of the subordinator, while the eigenfunctions remain the same. Thus the eigenvalues of the 
SubOU semigroup operator vf are {e -9 ^™)*, n = 0,1,...} with the same eigenfunctions. The 
general spectral representation of the semigroup of a subordinate symmetric Markov process 



can be found in Okura ( 2002 ) and in Alberverio and Riidiger ( 2003 ) and Alberverio and Riidiger 



(2005), where subordinate OU processes and their semigroups are studied in the general setting 



of symmetric Markov processes. Applications in finance can be found in Linetsky (2007), 



Mendoza et al. (2010) and Mendoza and Linetsky (2010). 



For t > the eigenfunction expansions on the RHS of (2.11) and (2.13) for the OU and 
the SubOU semigroup converge to Vtf and Vff in the L 2 -norm for any / G L 2 (M, m). In 
financial applications, we are interested in pointwise convergence, as we need to compute values 
at specific levels of the underlying variable. For t > pointwise convergence results are available 
for OU and SubOU semigroups. 

Theorem 2.7. (1) The eigenfunction expansion (2.11) converges to Vtf{x) pointwise in x for 
each t > and each f G L 2 (M,m). 

(2) If either of the following condition is satisfied: (i) f(x) = Yl^La fnfn(x) converges absolutely 
for all x G K, or (ii) YlnLi e~^ Kn ^ t n~ x l /!L < oo for all t > 0, then the eigenfunction expansion 
(2.13) converges to vf fix) pointwise for all x £i for each t > and each f G L 2 (M,m). 

The eigenfunction expansion (2.11) for the OU semigroup converges pointwise without any 
further conditions for each t > and / G L 2 (M, m). The eigenfunction (2.13) for the SubOU 
semigroup converges pointwise for each t > and / G L 2 (M,m) under the mild sufficient 
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condition on the Laplace exponent of the subordinator in (2) of Theorem 2.7. In practice 
this condition is satisfied for all subordinators with drift 7 > due to the factor e~~ /Kt . In 
the pure jump case 7 = 0, it is satisfied for all tempered stable subordinators with p > 0. 
Furthermore, for subordinators for which it is not satisfied, while the eigenfunction expansion 
(2.13) is not guaranteed to converge pointwise for each t > and each / £ L 2 (M, m), it may 
converge pointwise for some t > and some functions /, depending on the rate of decay of the 
coefficients f n as n increases. 

We also have the following expansions for OU and SubOU transition densities. 

Theorem 2.8. (1) The OU transition density (2.1) has the eigenfunction expansion 

00 

p(t, x, y) = m(y) £ e-^Vn^Kfo) (2.14) 

?1=0 

converging for all t > uniformly in x, y on compacts. 

(2) If the Laplace exponent of the subordinator satisfies Y^=i e~^ Kn ^ t n~^ < 00 for all t > 0, 
the SubOU transition density has the eigenfunction expansion 

00 

p*(t, x, y) = m(y) £ ^p n {x) Vn {y) (2.15) 

n=0 

converging for all t > uniformly in x, y on compacts. 

In the numerical implementation one needs to truncate eigenfunction expansions after a 
finite number of terms. Truncation error bounds of the expansion (2.15) in the I? and the 
pointwise sense can be easily derived. Here we present the pointwise error bound, as it is of 
most interest in finance. I? bounds can be derived similarly. 

Theorem 2.9. Suppose that the Laplace exponent of the subordinator satisfies J2^Lo e -< ^ Kn )* < 
00 for all t > 0. Then for any f £ L 2 (M, m), the truncation error has the following bound: 



e-^fnMx] 



n=M 



< 1.0864||/||e ^ e 



9 00 
(x-e 2 \ . 

<p(Kn)t 



n=M 



If 7 > 0, we can derive a particularly simple pointwise truncation error estimate: 



e-«" n *f n <p n {x) 



n=M 



< 1.0864 



1 - e~i Kt ' 



From these estimates it is clear that the convergence rate is governed by the OU mean reversion 
rate k and time to maturity t, as well as the Laplace exponent of the subordinator. The greater 
the k and the longer the time to maturity, the faster the convergence. In particular, if 7 > 0, the 
convergence is exponential. In the pure jump case 7 = with tempered stable subordinators 
with p > 0, the truncation error can similarly be shown to be 0( e *cr(- P )( K M)^ with r(-p) < 
and C > 0. We note that these error bounds are conservative since they rely on the estimate 
|/n| < 11/11- Depending on the properties of /, the coefficients /„ may converge to zero at a fast 
rate, resulting in faster convergence than is implied by these estimates. 
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3 Commodity Models With Mean- Reverting Jumps 

3.1 Futures Dynamics 

We start with (Q, F, (Ft)t>o) as in section 2.3 endowed with a probability measure Q and 
assume that, under Q, the canonical process X is a SubOU process with generating tuple 
(k, 8, a, 7, v) and starting point Xq = xo £ R. Let {F(0, t), t > 0} be the initial futures curve (a 
given deterministic function of time). We take Q to be the risk- neutral pricing measure chosen 
by the market and model the commodity spot price St under Q as the (scaled) exponential of 
the SubOU process X: 

S t = F(0,t)e Xt ~ G ® . (3.1) 

The function G(t) is selected so that the expectation of the spot price under Q is equal to the 
initial futures price, E Q [S t ] = F(0,t), which implies G(t) = \nE Q [e Xt }. 
To compute futures price dynamics, we need the following. 

Lemma 3.1. The expansion of the exponential function in the eigenfunction basis (2.13) reads: 

OO 2 1/ \ n 

e x = Y,fn<Pn(x), fn = e e+2 r~[^=) . (3.2) 



Jn\ Vv2k/ 
n=0 

The expansion converges absolutely for each 

We now compute the futures price process {F(s,t) = E®[St\ F s ], s G [0,t]} for each fixed 
maturity t > using Lemma |3.1| 



Theorem 3.1. (1) The function G(t) in the model (3.1) is given by: 

oo 

E Q [e x t] = J2e-' KKn)t fn<fn(xo), (3.3) 



•AO 

n=0 



where f n are given in (3.2), and the expansion converges absolutely for each xq £ R, all t > 
and any Laplace exponent (j). 

(2) For each fixed maturity time t > 0, the futures price F(s,t) is a martingale on [0,t] given 
by: 

oo 

F(s,t) = F(0,t)e- G ^^2e-^ n ^f n M^s), s G [0,t\. (3.4) 

n=0 



At time zero, s = 0, (3.4) reduces to the identity F(0,t) = F(0,t). At maturity, s = t, the 
futures price is equal to the spot price and (3.4) reduces to (3.1 ) due to Eq.(3.2). Eq.( |3.4[ ) gives a 
martingale expansion for the futures price. Note that for each n the process {e^ Kn ^ s tp n {X s ), s > 
0} is a martingale due to the eigenfunction property: 

E% n (X s )\X t ] = e-^"^-V(^)- 

Thus, Eq.(3.4) represents the futures price process as an expansion in martingales associated 
with the eigenfunctions of the SubOU semigroup. 

Since the process X can be expressed in terms of the spot price process S and the initial 
futures curve by inverting (3.1), 

X s = \n(S s /F(0,s))) + G(s), (3.5) 



Eq.(3.4) expresses the dynamics of the futures price in terms of the spot price dynamics and 
the initial futures curve. Alternatively, we can view Eq.(3.4) as the process for the futures price 
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driven by the SubOU process X without any reference to the spot price S. In this interpretation, 
our model can be viewed as the model for the evolution of the futures curve, rather than the 



spot price model. Eq.(3.4) directly defines the martingale futures dynamics. The spot dynamics 
(3.1) then follows as the limiting case. 

Remark 3.1. The Case without Time Change. When Xt is an OU rather than SubOU 
process, our model reduces to the standard exponential OU model: 

S t = F(0,t)e 



X t ~x e- Kt -e{l-e-^)-^-{l-e- 2 ^) 



By applying Ito's formula, we obtain the spot price SDE: dSt = «(©(£)— In St)Stdt+aStdB t with 
8(t) = I (| lnF(0, £) + g(l - e" 2K *)) + lnF(0,t). This is essentially the same SDE as the 
Schwartz ( 1997[ ) but with the long run level @(t) taken to be a deterministic function 



Model 1 in 



of time completely de termined by the initial futures curve. Using the generating function of 
Hermite polynomials ( |Lebedev| fl972| ) p.60), £n=o ^H n (z) = e 2zw ~ w \ when X t is an OU 



process (i.e., <fi(\) = A), Eq.(3.4) reduces to: 



( 2 

F(s, t) = F(0, t) exp \ Xse-^-^ - x e~ K * - #(e~ K( '- s) - e~ Kt ) - °— le' 2 "^ - e - 2Kt ) 

/ o \ exp{-«(t-s)} f 2 

= F( °' t) to) eXP {-4-^^' - ^ - 

This expression for the futures price dynamics in terms of the initial futures curve and the 
spot price dynamics in the OU model can be found in Clewlow and Strickland (1999), Eq.(2.5). 
Using Ito's formula, one can show that 

dF(s, t) = ae- K{t - s) F(s, t)dB„ s G [0, t]. (3.6) 

We now discuss futures dynamics under the physical measure P. The form for the futures 
process is still given by ( |3.4[ ). However, the law of X changes under an equivalent meas ure 

p p jr^ 

change. Let (B ,C ,11 ) be the semimartingale characteristics of X under P. Theorem 



2.5 



gives the general conditions on the semimartingale characteristics of (B ,C ,11 ). Any semi- 
martingale satisfying these conditions can be chosen as a candidate driver for the commodity 
model under P that leads to the model driven by the given SubOU process with generating tuple 
(k, 9, <T, 7, u) under Q. In order to retain analytical tractability under P, we are interested in 
equivalent measure transformations that transform a given SubOU process into another SubOU 
process plus possibly a deterministic function of time. Using Theorem 2.4 and Theorem 2.5, we 
obtain the following result. 

Theorem 3.2. Consider the canonical process X on (Q,F, (Ft)t>o)- Suppose under measure 
Q with Q(Xq = xq) = 1 the canonical process X is a SubOU process with generating tuple 
(k, 6, o~, 7, v) and under measure P with P(Xo = xq) = 1 it is a SubOU process with gener- 
ating tuple (kp, Op, o~p, 7p, vp) plus a deterministic function H(t). Then Q and P are locally 
equivalent if and only if: 

(1) H is absolutely continuous with H(0) = i/7 > 0, and H(t) = for all t j/7 = 0. 

(2) -f P ap = ja 2 . 

(3) the Bellinger condition f y ^ (\/ 7tP ( x , y) — y/ir(x, y)^ 2 dy < 00 is satisfied. 
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The Hellinger condition (3) can be simplified using Proposition 2.4 For example, in the 



case where the Levy measures v p and are both those of tempered stable subordinators, the 



Hellinger condition (3) reduces to the conditions presented in Corollary 2.2 

If X under F is specified to be a SubOU process plus some deterministic drift given by the 
function H(t), the model parameters can be estimated from the time series of futures prices 
by filtering methods. In this case the transition density of the underlying SubOU process is 



known explicitly and given by (2.15). The pure OU diffusion based model has been estimated 



by Schwartz (1997). In that case, the noise term is Gaussian and the standard Kalman filter 



can be used. In our SubOU case, the noise term for the transition equation is not Gaussian, 
and the particle filter algorithm (or the extended particle filter or the unscented particle filter) 



can be used since we know the transition density of X in closed form (see Haykin (2001) and 



Javaheri et al. (2003)). 



3.2 The Maturity Effect 



The maturity effect (also known as the Samuelson hypothesis, see Samuelson (1965)) in the 
commodities futures markets is the well-known increase in commodity futures price volatility 
as the futures contract approaches maturity. The maturity effect implies that long term futures 
are less volatile than short term futures, and is well documented in the empirical literature (see 
Bessembinder et al. (1995), Kalev and Duong (2008) and references therein). The maturity 
effect is obviously present in the pure OU model (3.6), where futures volatility ae~ KT decays 



exponentially as time to maturity r = t — s increases, with mean reversion rate controlling the 
rate of decay. Here we investigate the maturity effect in our SubOU model. 

We start with characterizing futures volatility in the general semimartingale setting. For 



a futures contract with maturity time t, define 



111 F(Q,t) ■ 



s £ [Q,t], the cumulative contin- 



uously compounded return process over the time interval [0, s] with s < t. Since F(s, t) is a 
semimartingale, r* is also a semimartingale. We measure volatility of the futures return process 



experienced over the time interval [0, s] by its quadratic variation (QV) 



r , r 



(the square- 



bracket process). This definition of volatility has been widely used in the econometric literature 



(see Andersen et al. (2009)). With this definition, the maturity effect can be mathematically 



defined as follows. 

Definition 3.1. A futures model is said to exhibit the maturity effect almost surely if 

P([r* 1 ,r* 1 ] s > [r t2 ,r t2 ] s ) = 1 for any < s < h < t 2 . 
Remark 3.2. If IP and Q are locally equivalent, then the QV of a semimartingale under P is 



a version of the QV under Q (Jacod and Shiryaev (2003) Theorem III. 3. 13). Hence, if the 



maturity effect is present in the futures dynamics under the physical measure, it is also present 
under the risk-neutral measure. We will compute [r , r*] s under Q. 



Remark 3.3. In the pure OU model (|3.6|) the QV of futures return process is a deterministic 

^2 



function [r*, r*] { 
Note that [r*. r* 



-2nt ( „2ks 



2k 



1) decreasing in t for each fixed s, < s < t. 



„tc ~icl 



] s + ^ u<s (Ar*) , where r denotes the continuous martingale 



part of the process r t . From Eq.(3.4), r* = —G(t) + 8 + ? — h lng(X s , s, t), where the function 



9 is: 



g(x,s,t) := J2 e ~ 



(ren)(t— s) 



n=0 



a 



—.tl r , 



a 



4 



Since we know the semimartingale characteristics of the SubOU process X, from Kallsen (2006) 
Proposition 2.5 we know that [r tc ,r% = 7 a 2 Jo f ^%^ ) : ''du = 7a 2 / p s ( 9 *£ n -?$ ) du 



1 g(X u .,u,t) 
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and (Ar*) 2 = (r* -r*_) 2 = (]ng(X u ,u,t)-lng(X u _,u- ,i)) 2 = fgf ds)'. There- 

fore, [r*. r% = 7 , 2 £ (fg^f ) 'du + E u < s ( J™ 1 : 2 . Note that ,(*, «, t) > 

and g x (x, u, t) > for each x, and g(x, u, t) depends only on t — u. We thus have the following 
result. 

Theorem 3.3. If g x (x, 0, t)/g(x, 0, t) is decreasing in t for each x, then the maturity effect holds 
in the SubOU model. 

While the condition in Theorem 3.3 is hard to check analytically since the function g is 
given by the Hermite expansion, it can be easily checked numerically. We carried out extensive 
numerical testing for a wide range of parameter scenarios in pure jump (7 = 0) and jump- 
diffusion (7 > 0) cases and verified that it was indeed satisfied in all the cases. We thus 
conjecture that the condition in Theorem 3.3 is satisfied, and the maturity effect holds for our 
SubOU models. 

Figure [2] illustrates the maturity effect as follows. We simulated 10,000 sample paths on the 
time interval [0, 1/2] of pure jump (7 = 0) SubOU processes X with parameters 6 = 0, a = 0.5, 
with the Inverse Gaussian subordinator with mean rate \x = 1 and variance rate v = 1, and with 
k = 0.01,0.1, and 1. We then constructed 10,000 sample paths of the futures price processes 
with maturities 1/2, 1, 2, 3, 4 and 5 years for each of the underlying SubOU processes using the 
model relationship ( |3.4[ ) under Q, estimated realized quadratic variations of futures returns on 
each sample path (the quadratic variation is the same under P and under Q), and verified that 
[r* 1 , r t:L ]o.5 > [t* 2 , t* 2 ]o.5 for t\ < t<i on each sample path. Figure [2] plots the estimated mean of 
the quadratic variation of futures returns as functions of futures contract maturity for the three 
values of the rate of mean reversion n = 0.01,0.1, and 1. The maturity effect is clearly seen in 
the plot. As in the pure diffusion OU model, k controls the maturity effect in pure jump and 
jump-diffusion SubOU models. 




Figure 2: E^{[r*, r*]o.s} as a function of futures maturity t for a SubOU Process with Inverse 
Gaussian subordinator (with parameters k = 0.01,0.1,1, 9 = 0, a = 0.5, 7 = 0, mean rate 
(j, = 1, variance rate v = 1). 

To further illustrate, Figure [3] plots a sample path of the driving SubOU process in the jump- 
diffusion case and the corresponding futures price process with 3 years to maturity at time zero. 
The maturity effect is clearly seen in the sample path dynamics, as the futures price experiences 
low realized volatility far away from maturity, and the realized volatility substantially increases 
as the futures contract approaches maturity. 
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(a) the SubOU process 



(b) the futures process 



Figure 3: A Sample Path of a SubOU Process and the corresponding futures price process with 
the Inverse Gaussian subordinator (with parameters k = 1, 9 = 0, a = 0.5, 7 = 0.1, mean rate 
[i = l, variance rate v = 1). 



Remark 3.4. It is important to note that the rate of mean reversion k that enters the expression 



for the diffusion volatility in the pure diffusion OU case (3.6) and in the quadratic variation 



process through the functional form (3.4) of the dependence of the futures price on the SubOU 
process in the SubOU model is the rate of mean reversion under the risk-neutral pricing measure 
Q. It is the risk-neutral rate of mean reversion that controls the maturity effect. That is, the 
presence of the maturity effect in the futures time series under the physical measure P is governed 
by the rate of mean reversion under the pricing measure. If there is no mean reversion under 
the pricing measure Q, i.e., X is taken to be a subordinate Brownian motion under Q rather 
than a subordinate OU process, there is no maturity effect under P. Thus, the presence of the 
maturity effect under P requires X to be a SubOU process under Q, as futures models built on 
subordinate Brownian motions (Levy processes) do not possess the maturity effect. In contrast, 
SubOU models are capable of modeling the maturity effect. 

3.3 Futures Options Pricing 

We consider pricing European put and call options on a futures contract. Suppose the strike 
price is K. The underlying futures contract matures at time t* and the option expires at t < t*. 
The time r = t* —t varies across commodities, ranging from several days for natural gas to one 
month for gold. 

Here we only consider pricing the put option. The call option price is given by the put-call 
parity. Alternatively, a similar eigenfunction expansion can be obtained for the call pricing 
function, and the put-call parity can be verified directly. The put payoff at expiration t is 
(K — F(t,t*)) + , where F(t,t*) is the immaturity futures price at time t. In our model it is 



related to Xt by (3.4). It is convenient to write the payoff function as follows: 

(K - F(x, t, t*))+ = {K- F(x, t, t*))l {x<x * } , 

where x* is the unique solution of the equation F(x, t,t*) = K, and F(x, t, t*) is the immaturity 
futures price at time t as a function of the state variable Xt = x given by ( |3.4[ ). Since F{x, t, t*) 
is a strictly increasing function of x, the solution to this equation is unique and can be easily 
computed numerically using bisection or any other root bracketing algorithm. To price the put 
option at time zero, we thus need to first find x* corresponding to the strike price K and then 
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compute the expectation in: 

P(t, t*,F(0, t*),K) = B(0, t)E® [(K - F(X t , t, t*))l {Xt<x * } 

where B(0, t) is the risk-free discount factor from the option expiration t to time zero. 

Theorem 3.4. Let x* be the unique solution of the equation F(x,t,t*) = K and define w* := 
^(x* — 0), t := t* — t, a := an d F := F(0,t*). Suppose the Laplace exponent of the 

subordinator satisfies Y2n=i e -< ^ Kn )*n~3 < oo. Then the put price has the absolutely convergent 
eigenfunction expansion: 



P(t, f,K, F) = B(0, t) e-^pnit, t* ,w* , F)ip n (x ), 



n=0 



(t, t*,w*,F) = —^— Ubniw*) - Fe* 
V7r2 n n! 



2 00 

fs-G(f ) e 

m=0 



-<f>(Km) 



.a' 



ml 



(3.7) 



.(»*)}■ (3-8) 



b n (w) 



H n (x)e x dx 



H n -i{w)e~ 

min(n,m) 



n = 0, 

n = 1,2,--- 



H m {x)H n (x)e- x2 dx= (TJIJ 2 ^ 6 
■°° k=o V / V / 



'n+m— 2k 



(w). 



(3.9) 



(3.10) 



The call price is given by the put-call parity C(t, t*,K, F) = B(0, t)(F - K) + P(t, t*,K, F). 



Remark 3.5. The option written on the spot price is obtained by setting t = t* in (3.8). 



Remark 3.6. The Case Without Time Change. In the pure diffusion OU model, the option 
pricing formulas collapse to the Black-Scholes type formulas for the exponential OU diffusion 



model obtained by Clelow and Strickland (1999): 



P(t,t*,K,F) = B(0,t) [K$(-d-) - F<S>(-d + )\ , C(t,t*,K,F) = 5(0, t) [F$(d+) - K$(d-)] , 



M£)-ge- 2Kr (i-e- 



■2nt\ 



'2k 



,-2/tt 



4 Stochastic Volatility and Time Inhomogeneity 

Models based on SubOU processes described in the previous section can be calibrated to 
fit a variety of volatility smile patterns observed in commodity options markets. However, 
they are generally not flexible enough in order to fit the entire volatility surface across different 
maturities. In this section we study a further extension of SubOU models to introduce stochastic 
volatility and time inhomogeneity, such as seasonality in options' implied volatility typical for 
some commodities, such as natural gas. 

We consider absolutely continuous time changes of the form 

T t = J (a(u) + Z u y u , (4.1) 

where a(t) > is a deterministic function of time and Z is a CIR diffusion solving the SDE 

dZ t = k z (6 z - Z t )dt + az\fZ t dB u Z = z , 
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with parameters assumed to satisfy the Feller condition, 2QzKzjo\ > 1 to ensure that zero is 
an inaccessible boundary. 

The activity rate process a(t) + Zt has the form of the so-called CIR++ process well known 



in the interest rate modeling literature (e.g., Brigo and Mercurio (2006)). The advantage of 
the CIR process is in its analytical tractability. Its transition probability density, the Laplace 
transform of its integral, and the Laplace transform conditional on the terminal state of the 
process are all known in closed form. The relevant results are collected in Appendix |A} 

Define the process S to be the inverse of T, St ■= inf{u > : T u > t}. Since T is a strictly 
increasing continuous process, so is S. It is also clear that T t = inf{u > : S u > t}. 

Assume that on some complete probability space (f2, J 7 , P) we have a cadlag SubOU process 
X with generating tuple (k, 9, a, 7, v), Xq = xq and an independent absolutely continuous time 
change T of the form in (4.1). Let (J r t)t>o be the smallest right-continuous complete filtration 
generated by the processes Xt, Z$ t and St- Then Tt is a stopping time w.r.t. {Ft)t>o for every t, 
and we can define the time changed filtration Qt ■= J-T t - R is clear that T and Z are adapted to 
{Qt)t>Q- Define a new process Y by Yt := Xx t , with lb = 2/0 = xq. From Jacod (1979) Corollary 
10.12, Y is a (^)-semimartingale, and from Kallsen and Shiryaev (2002) Lemma 5 it admits 
the following local characteristics (B,C,H): 



t r roo 

(o(«) + Z,(w)) 7 K(0-y s _(a;))+ / 
L Jo 



\x\<l 



xp(u;Y s -(uj),Y s -(lj) + x)dxv{du) 



ds, 



Ct(oj) = 7<j / (a(s) + Z s (ui))ds, II(cj, dt, dx) = (a(t) + Zt(uj))ir(Yt-(u), x)dx, 
Jo 



where 7r(-, •) is defined in Theorem 2.1 and here x is interpreted as the jump size. From these 
expressions we see that the role of the absolutely continuous time change is to scale all the local 
characteristics of the SubOU process with the stochastic activity rate or stochastic volatility. 
The bivariate process (1", Z) is also a (^)-semimartingale. We have the following result on its 
cross- variation process. 



Proposition 4.1. The cross-variation process [Y c ,Z c ]t 
ous local martingale parts of Y and Z respectively. 



0, where Y c and Z c are the continu- 



It is clear that (Y,Z) is also a Markov process w.r.t. the filtration {Qt)t>o- Given Qt, the 
distribution of Y t+S depends only on T i+S — T t and Yt, and T t+S — T t depends only on Z t . The 
distribution of Zt+ S depends only on Zt. Thus, conditional expectations of the form E[/(Yt)|£7 s ] 
reduce to E[/(F t )|Y s , Z s ] by the Markov property. Using conditioning and the spectral represen- 
tation of the SubOU semigroup, such expectations can be computed in terms of eigenfunction 
expansions. 



Theorem 4.1. For f G 1? 



m), suppose one of the following two conditions is satisfied: 

Y^=ofn^Pn{x), where f n = (f,(p n ), converges abso- 



(1) The eigenfunction expansion f(x] 
lutely for each x. 

( 2 ) T.n=o e~^ Kn ^s a{ ^ du Ccm (t - s, <f>(Kn) \z)n~i < 00 for some z>0 (and hence for all z; 
it is straightforward to show this using (A.4) ), where the Laplace transform Cqir {t, ■ \z) 
is given in Appendix A. 



Then E[f{Y t )\Y s ,Z s ) = ££L e"^) £ a ^ du £ C iR (t - s,^n) Z s ) f n ip n (Y s ) . 
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We can now introduce stochastic volatility and time inhomogeneity in commodity models. 
Let Yt = Xx t be the time changed SubOU process as above. Under the risk-neutral pricing 
measure Q chosen by the market, we model the spot price as follows: 



Y t -G(t) 



S t = F(0,t)e 

where the function G(t) is selected so that e G ^ = E 1 ® 
exponential function, we obtain the futures price process. 



on- 



Applying Theorem 



4.1 



(4.2) 
to the 



Theorem 4.2. (1) e 



G(t) 



Y,n= z-* {Kn)S ° a{u)du £cm(t, 



(an) z jf n ip n (y ), where f n are 
given in Lemma \3.1\ The expansion converges absolutely for all zq > 0, yo G M, and any 
Laplace exponent <fi. 

(2) For each t > 0, the futures price is a martingale on [0,t] given by: 



F(s,t) = F(0,t)e- G ^£< 



-4>(ku) f 



n=0 



a{u)du C C m(t 



Z s )f nVn (Y s ), 8€[0,t]. (4.3) 



To investigate the maturity effect, we need to compute the QV process [r*,r ] s , which is 



more involved in this case due to the extra state variable Z. From (4.3), r* 
lng(Y s ,Z s ,s,t), where 



g(y, z,s,t) = f2 e~ HKn) J * a(u)du C C iR (t - s, cp{nn) 



n=0 



a 



n 1 



(y-o) 



Again we use Kallsen (2006) Proposition 2.5 to compute [r ,r ] s from the local characteristics 
of the semimartingale (Y,Z). Since the cross- variation is zero by Proposition 4.1 we do not 
have cross derivative terms and obtain: 



r , r 



[r tc ,r% +1 a 2 



(o(«) 



+ 



+*)( »ff-'f"*' ) )'*.+.* fz. 

\ g(Y u _,Z u ,u,t) ) J 
YuVYu - g y (y,z,u,t) 



g z (Y u _,Z u ,u,t) \' 
g(Y u -,Z u ,u,t) ) 



du 



E 

u<s 



y u ay u „ g(y,z,u,t) 



dy 



Note that g and g y are positive, but g z is not necessarily so, and g(y, z, u, t) depends on u and 
t only through t — u. It is thus clear that we have the following: 

Theorem 4.3. If ~7~~g^r and f^^ny) o r e decreasing in t for any (y, z), then the maturity 



effect holds. 



g(y,z,o,t) 



As in the SubOU case in section 3.3, this condition is hard to check analytically, but can 
be easily verified numerically. We have conducted extensive numerical experiments and verified 
this condition for all parameter specifications we have tested. 

For the model with stochastic volatility the option pricing formula is more involved since 
the futures price F(t,t*) at expiration of the option t is now determined by the values of two 
state variables Yt and Z± at that time, F(t,t*) = F(Yt, Zt,t,t*). We condition on the state of 
the CIR process Zt at time t and reduce the problem to the SubOU case. One then has to 
use the conditional Laplace transform (A. 5) instead of (A. 3), since we have conditioned on Z%. 
Hence, the pricing formula is expressed as an integral with respect to the transition density of 
the CIR process (A. 2). An additional subtlety is that y* now depends on z. Namely, for each 

y*{z) such that F(y*,z,t,t*) = K. Then the put payoff 
(K - F(y,z,t,t*))l {y<r{z)} . 



fixed z > 0, there exists a unique y 
function can be rewritten as (K — F(y, z, t, t*)) 
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4.1 



Theorem 4.4. For each fixed z > 0, let y*{z) denote the unique solution of the equation 
F(y, z,t,t*) = K, where F(y, z,t,t*) is the futures pricing function (4.3). Define w*(z)^ 
^(y*(z) — 9), t := t* — t, a := and F := F(0,t*). Suppose condition (2) of Theorem 
and the following condition are satisfied: 

oo 

J2e~^ Kn) ti a ( u ) du £ CIR (t,<p(\) \zo,z)n-* < oo for some z and hence for all z. (4.4) 

n=0 



(It is easy to show this using (A.6).j Then the put price is given by: 

P(t,t*,K,F) = B(0,t) 



(4.5) 



x / i ]>>~ 0(ren)/ ° a(u)du ^^ 



v n=0 



where pcm(t, zq, z t ) is the CIR transition density (A. 2) and 

Pn (t,t*,W*{z),F) 



1 



Vvr2 n n! 



(4.6) 



2 t* n m 

x { Kb n (w*(z)) - Fe e+ ^- G ^ V e ~^ m )/* a ^ du C C m(r, <j>( Km )\z t )— a n , m (w*(z)) 



m=0 



where b n (w) and a n ^ m {w) are given by (3.9) and (3.10). The call price is given by the put-call 
parity. 

Remark 4.1. For options written on the spot price, in contrast to futures options, we only 
need the Laplace transform of the time change instead of the conditional Laplace transform. 

K 



Furthermore, in this case y* = In y F ^ + G(t) is independent of Z t . By setting t = t* and 
using C C m(t, (j)(Kn)\z , zt)pcm(t, z , z t )dz t = C C m(t, cp(Kn)\z ), the put price becomes 

oo 

P(t,K,F) = ^ e -^")/o«H^£ c/R (t,0( K n)|z o K(t,t,^,F)^n(yo). 

n=0 

5 Model Implementation and Calibration Examples 

The models introduced in this paper were implemented in C++ on a PC. Hermite expansions 
can be efficiently computed using the following classical recursion for Hermite polynomials 
( Lebedev| ( [19721 ) p.61): 

H (x) = l, H 1 (x) = 2x, H n (x) = 2xH n _ 1 (x)-2(n-l)H n ^ 2 {x), n>2. 



To compute the option pricing formula (3.7), we need to evaluate the coefficients b n and a n , m . 
From the equation (3.9), it is easy to see b n can be computed recursively using the recursion 
for Hermite polynomials. Equation (3.10) is a closed-form formula for a nm , but it is not 
convenient to use from a computational perspective. We have the following computationally 
efficient approach for a n ^ m . 

Proposition 5.1. The coefficients a njTn satisfy the following: 

ao,o( x ) = V^^(V2x), a n>n (x) = 2na„_i ;n _i(x) - H n - 1 (x)H n (x)e~ x2 , n > 1, (5.1) 

a nm (x) = Hn{x)H m+1 (x) - H (x)H n+l {x) e _ x ^ n ^ m> n > q, m > 0. (5.2) 

2(m — n) 
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To evaluate the option pricing formula (4.5) for the model with stochastic volatility, we first 
truncate the integral in zt at some level M large enough that the probability of the CIR process 
to exceed M at time t is less than the desired error tolerance. We then use the Simpson rule to 
discretize the integral on the interval [0, M\. The CIR transition density at each node zt{k) is 
computed by (A. 2), at each integration node zt(k) the value of y*{k) is found by the bisection 
algorithm, and the integrand is computed similar to the option pricing formula (4.5) in the 
SubOU case (with the distinction that under the time changed SubOU the conditional Laplace 
transform (A. 5) enters the expression in place of the Laplace transform (A. 3) in the SubOU 
case) . 

CPU times generally depend on time to maturity t and the model parameters. For short 
maturities (say, less than two weeks to expiration), one may have to use infinite-precision 
arithmetics to achieve required accuracy in summing up the series. To compute short maturity 
option prices we used the GNU MP Bignum library. For longer maturity options double precision 
is sufficient. In our numerical experiments on a PC running Linux (Intel Core 2 Duo CPU at 
2.53GHz with 2.00GB RAM), CPU times ranged from several milliseconds up to hundreds of 
milliseconds per option for the SubOU model, depending on the combination of parameters, 
and from hundreds of milliseconds up to several seconds per option for the SubOU model with 
stochastic volatility. 

We now present calibration examples of the SubOU model with the IG subordinator to 
implied volatility smile curves extracted from market prices of options on six commodity futures. 
We have also calibrated for other commodities, and the results are similar to what are displayed 
here. However, due to space constraints, only six of them are shown. The commodities included 
two metals (copper, gold), two energies (crude oil, natural gas), and two agricultural (corn and 
wheat). Market data on implied volatilities for this study were provided by Morgan Stanley's 
Commodity Stategies Group and were extracted from commodity futures options market prices 
on July 2nd 2009. All options had approximately six months to expiration. The moneyness 
defined as the ratio of the option strike price to the futures price ranged from 0.6 to 1.8 for all 
commodities. To calibrate the model to market implied volatilities, we minimized the sum of 
squared differences between the market and the model implied volatilities. There are a total 
of six parameters in the SubOU model: three parameters of the background OU process and 
three parameters of the inverse Gaussian subordinator with drift. Without loss of generality, 
the starting SubOU state xq was set to zero (it can always be set to zero by changing 9 to 
9 — xq without affecting the option price) . Our calibration results are presented in Figure [4j In 
these instances the SubOU model with the IG subordinator provides an excellent fit to volatility 
smiles for all eight commodities (well within the bid/ask spread for each option). 

While the SubOU model calibrates well to commodity volatility smiles for a fixed matu- 
rity, it may generally lack flexibility to capture the entire volatility surface across both the 
maturity dimension and the strike (moneyness) dimension. The time changed SubOU model 
with stochastic volatility and possible time inhomogeneity has additional flexibility to capture 
time dependence in the shape and steepness of the volatility smile and time dependence in the 
at-the-money volatility term structure. In Figure [5] we calibrate the SubOU model with the 
inverse Gaussian subordinator time changed with the integral of the CIR process to the implied 
volatility surfaces for zinc. We used four maturities (6 months, 1 year, 1.5 years and 2 years) 
in our calibration. The deterministic activity rate component was taken to be a piecewise con- 
stant function (constant between adjacent futures maturity dates). The time changed SubOU 
model provided an excellent fit to this volatility surface (well within the bid/ask spreads for 
all options). The deterministic activity rate allowed us to capture the sharp decay in the ATM 
implied volatilities, the IG subordinator allowed us to capture steep smiles for shorter-dated 
maturities, and the CIR stochastic volatility supported the longer-dated smiles. In contrast, 
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ure 4: SubOU Model Calibration Results to Implied Volatility Smiles for Commodities 
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SubOU models without stochastic volatility exhibit faster flattening of the volatility smile as 
we go further out in maturity. 



Zinc 

Imarket 
[ ;model 




Figure 5: Time Changed SubOU Volatility Surface Calibration Results For Zinc 



6 Conclusion 

This paper studied a class of subordinate OU processes, their sample path properties, equiv- 
alent measure transformations, and the spectral representation of their transition semigroup. 
As an application, we constructed a new class of commodity models with mean-reverting jumps 
based on subordinate OU process. Further time changing by the integral of a CIR process plus 
a deterministic function of time, we induced stochastic volatility and time inhomogeneity in the 
models. We obtained analytical solutions for commodity futures options in terms of Hermite 
expansions and showed that the models exhibit the maturity effect and are flexible enough 
to capture a wide variety of implied volatility smile patterns observed in energy, metals, and 
agricultural commodities futures options. 

We are currently developing computational methods for American-style futures options in 
these models. It turns out that the eigenfunction expansion approach to pricing European 
options followed in this paper can be extended to Bermudan-style options with a finite number 
of exercise opportunities. Richardson extrapolation can then be used to obtain solutions for 
American-style options. 

In future work we plan to extend this class of models to multi-commodity products, such 
as spread options, and to path-dependent options such as Asian-style options. An extension 



to American-style options is developed in Li and Linetsky (2011). We also anticipate that 



subordinate OU processes studied in this paper will find other applications beyond commodities, 
such as in interest rate modeling, volatility modeling, and real options. 
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A CIR Processes 

Let {Zf, t > 0} be a CIR diffusion starting from Zq = z > and solving the SDE 

dZ t = k(6 - Z t )dt + ay%dB t . (A.l) 

Assume the long run level 9, the rate of mean reversion k, and the volatility parameter a satisfy 
the Feller condition d := > 1 to ensure that the process stays strictly positive (zero is an 
inaccessible boundary) . 

The CIR transition density pcm(t, zq, z) is given by 



Pcm(t,z ,z) 



2 k 



o- 2 (l - e~ Kt ) 



CT 2 (l-e- ret ) 



z e 



— Kt 



d-1 
2 



\k\[zq 



ze 



— Kt 



a 2 (l - e~ Kt ) 



(A.2) 



where Id-i(-) is the modified Bessel function of the first kind of order d — 1. 



The Laplace transform C C m(t, X\z ) := E Zo e~ A io z ^ du 
formula for the short rate process XZt. 

CciR(t,X\zo) = C(t,X)e- B ^ zo , 
where C(t,X) = — — J \/ .... — , B(t,X) = 



is given by the CIR bond pricing 

(A.3) 

2A(e^ A ) i - 1) 



( 7 (A) + K)( e 7Wt-l) + 2 7 (A)y '~^~> ( 7 ( A ) + K )( e7 (A)*_i) + 2 7 (A)' 

and 7(A) = V k 2 + 2a 2 X. The function Ccm(t, X\zq) has the following asymptotic behavior as 
A —7- 00: 

Ccm(tA\zo)-e W i-^V2Xt-^V2x\. (A.4) 



E 



The Laplace transform conditional on the state of th e process at time t, Ccni(t, X\zq, Zf) : 
e -\f*Z u du \z t = z t 



is also known in closed form (Broadie and Kaya 



(2006)): 



£cm(t,X\zo,zt) 



7 ( A ) e -0.5( 7 (A)-«)t^ _ e -,rf) 



k(1 - e-TW) 



x exp 



z + z t I k(1 + e" 



7 (A)(l + e- 



- 7 (A)t> 



Id- 



A~f(\)^fz^zi e -0.5 7 (A)t 
?2 1 _ e - 7 (A)t 



l-e 



-Kt 



1 _ e --r(x)t 



4K v / 5()2t e ~0.5Kt 



a 2 1— e Kt 

The function CciR{t, X\zq, zt) has the following asymptotic behavior as A — > 00: 
£ C /i?(t,A|z ,^) ~ Aiexp {-— V2Xt - ^±^\/2A 

B Proofs 



(A.5) 



(A.6) 



Part (2) of Theorem\2J^ Denote the RHS of Q*> in Theorem [2I] by Q* . If X' admits 
characteristics (B',C, II'), then from Ito's Formula for semimartingales, for any / G C 2 (M.) 



M t := f(X' t ) - f(x) - f G*f{X' s _)ds 
Jo 



is a local martingale. Since G^f £ Co(M) (the space of continuous functions on R vanishing 
at infinity), Q&f is bounded. f(X' t ) — f(x) is also bounded for all t. Hence E[M^] < 00 
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(M£ := sup s<t \M S \) for all t, and M is a martingale by |Protter (2005) Chapter 1 Theorem 
51. Note that from Theorem^ C£(R) 

is a core of D{Q^). Hence applying [Ethier and Kurtz 



(1986) Chapter 4 Theorem 4.1 to the martingale problem C 2 (IR)),P X ) and Corollary 4.3, 



it follows that 
Theorem 



x'- 



on the Skorohod space (fi, J 70 ). 



□ 



2.3 



If a; > 0, then for any y > 0, | -y + (x-9)(l-e~ Kt )\ < \y+(x—9)(l-e 
From (2.1 ), this implies p(t, x, x — y) > p(t, x, x + y) for any t > 0. Hence from the definition of 
7r(x, •), 7t(;e, — y) > 7r(x, y) for any y > 0. By integrating ir(x, •) on (— oo, — y) and (y, oo), we also 
get H(x, (— oo, — y)) > n(x, (y, oo)). The cases with x < 6 and x = 8 are proved similarly. □ 
Theorem \2-4[ Sufficiency. By replacing the Levy measure used in Remark 33.3 of 
Sato| ( |1999 ) by our state-dependent Levy measure, we can show that the Hellinger condition 



ly^o {V 7r '( x ^y) - V^i^y)) dy <oo implies that 



J\y\<i 



|y| • |7r'(X s _(a;),y) - ir(X s -(u), y)\dy < oo, 



(B.l) 



so h(x)(Y — 1) * n is finite. We first show that 



yp\u]X 8 -(u),X 8 -(u) +y)dyi/(du)- / yp{u;X s -{uj),X s -{ui) + y)dyv{du) 

[0,oo)J|y|<l J[0,<x,)J\y\<l 

Tr'(X s -(u),y) - ir(X s -(u),y) dy. 



y 

'\y\<l 
Note that 




lim 

n— >oo 



yp'(u, x, x + y)dyv'(du) 

'[0,oo) -/|i/|<l 

yp'(u, x, x + y)dyu'(du 

[0,oo) ^l/n<[j/|<l 



yp(ti, x, x + y)dyu(du) 

yp(u, x, x + y)dyu(du) 



l/n<\y\<\ 



M<1 



[0,oo) Jl/n<|»|<l 
(vr'(x,y) - ir(x,y))dy, 



where the last equality comes from the Dominated Convergence Theorem since we have (B.l). 
So for all o> we have B' = B + 7a 2 /3 • t + 7t(x)(Y - 1) * n, C" = C, n' = Y ■ U. 



Since n(w, t, dx) = 0, it is clear that ajs = 00, where ajs is defined in Jacod and Shiryaev 



(2003) (JS) III. 5. 6. The process H defined in JS III. 5. 7 becomes the following in our case: 



H t (u) 



\a 2 (3 s (co) 2 ds+ f [ (V7r'pf s _(a;),y) - yJir(X a -{u),y))' ''dydi 

JO Jy^O V 7 



It is clear that the integrand in the above expression is cadlag for every u. This implies that 
Ht{uj) < 00 for every uj and t. Hence the process H does not jump to infinity as defined in JS 
III. 5. 8. This fact together with ajs = 00 implies that Hypothesis III. 5. 29 of JS holds. 

P' is the unique solution to the martingale problem (a(Xo), X\F' Q , B' ,C ,11'). As remarked 
before, local uniqueness also holds. Note that P =<! Po- Now all conditions stated in JS Theorem 
III. 5. 34 are satisfied, which implies P' =3! P locally. By interchanging the role of P' and P and 
similarly defining f3' and H' , we can prove Po =3! P implies P' =3! P locally. Hence P' ~ P locally. 

Necessity. If P' ~ P locally, then (i) holds, and (ii) is implied by JS Theorem III. 3. 24. 
The uniqueness of the solution to the martingale problem (ct(Xo), X|Po, B, C, n) implies the P- 
martingale representation property w.r.t. X (JS Theorem III. 4. 29), hence JS Theorem III. 5. 19 
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holds, which further implies JS IV. 3. 32. Now the conditions in JS Theorem IV. 3. 35 are satisfied, 
and this theorem implies that the Hellinger process of order \ (see JS Definition IV. 1.24) is 
given by 



1 



7 (j 2 /3 s (w) 2 ds + - 



Jy^O 



^y/n'(X s -(uj),y) - y/Tv(X s -(u}),y)^ dyds 



JS Theorem IV. 2.1 says that h^{u) < oo both P and P'-a.s., hence there exists xq G M. such that 
fy^o (^\/ ^'(^O) y) — \/n(xo, yjj dy < oo. But one can show that the tail behavior at y = of 
\J tt'(x, y) — \/ir (x, y) does not depend on x (see Proposition 2.3, whose proof does not depend 

on Theorem 2.4), so we have J y _^ (^Jtt'(x, y) — \J ir(x, yfj dy < oo for any x. 

Therefore, the process H defined in the proof of the sufficiency part does not jump to infinity. 
This together with ajs = oo allows us to apply JS Corollary III. 5. 22 (ii) which gives the form 
of the density pr ocess . □ 

"'" '' — | , the transition density of 



Proposition 



2.2 



Define q(s, 0, y) :- 



V TT(T 2 , 



exp 



2<T 2 S 



Brownian motion starting at with drift k(6 — x) and volatility a. It is easy to see that 

p(s,x,x + y) = 
s^o q(s,0,y) 

uniformly for y on any compact interval. We wish to prove that 

J[ ,oo)P( s > x > x + y) u ( ds ) _ 



lim 



Note that 



J[o,oo) i( s > °' y) u ( ds ) 

f [os) p(s,x,x + y)v(ds) 



lim 



J[o,oo) P( s ' x -> x + y) u ( ds ) 



(B.2) 



(B.3) 



for any 5 > 0. This is because for s > 5, p(s, x,x + y) is bounded in s, and f< s ^ v[ds) < oo, so 
applying the Dominated Convergence Theorem 



lim 

y^O 



[<5,oo) 



p(s, x, x + y)v(ds) 



lim p(s, x, x + y)u(ds) 
[S,oo) 

exp{ 



(8-x) 2 (l 



which is finite, and hence lim,, 



|^(o?i 



y ~*° I[o, a o)P( s ' x ' x +y) L '( ds ) 

is replaced by q(s, 0, y) for the same reason. 



0. (B.3) is also true when p(s,x,x + y) 



Fix an interval [-M, M] for y. Then for any e > 0, there exists some 5 > 0, such that for 

(s,x,x+y 
q(s,0,y) 



any y £ [-M, M], 1 — e < p ^t^ <l + eifs<5. Hence 



for any y 6 [— M, M]. Now letting y — > we have 

z r I m P(^ x ^ + y)v{ds) 
1 — e < lim — < 1 + e. 
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. . J*,„ s , q(s,0,y)v(ds) 

Equation ( |B.3[ ) and lim^ ^ ' ^ g(sAy)l , ((fa) = 1 imply that 



lim — — — = lim — ' 

/[o )00 ) °' yM^s) J[ ,<5)9('5,o,y)^(ds) 

/[0 jO o)*'( s . !E > aj + l'M d6 



Hence 



1 — e < lim 



J[o,oc) ^(s, 0, 



< 1 + e 



for any e. Now letting e — > 0, (B.2) is proved. 



□ 



Proposition 2.3. Now we prove that the asymptotic of the Levy density 7f(y) of a SubBM 
does not depend on the drift. Suppose the drift and diffusion coefficients are n and a respectively. 
Then we have 

1 r (y - us) 2 



V2TTCT 2 : 



exp 



•}i/(ds). 



Similar to the proof in Proposition 2.2, it is straightforward to show that 



lim ir(y) = lim / 
which does not depend on \x. 



1 



[0,oo) TTO^S 



exp 



2a 2 , 



-\v{ds), 

• J 



□ 

Proposition 2.4 We prove the case with condition (1) here. The case with condition (2) 



is proved in |Kim et al. (2010). We can write: 



?r(y) 



1 

2cr 2 £ 



V2na 2 , 



k- y/2lt(J 2 S 



exp {-2k} Ks)ds - 



Similar to the proof in Proposition 2.2 the second integral on the RHS is finite as y — > 0, so we 



only need to be concerned with the first integral. The rest of the proof is similar to Song and 
Vondracek| (|2009[). Let u = y 2 /(2a 2 s). Then 



i 

2<T 2 { 



1 

V2na 2 ; 



exp 



{ - ^h s)d& 



\y\ 



2a 2 ^ 
{2a 2 f~ 1 



°° 3 V 2 



, ..'''•'■in-,, 



U 2 g 



where h(y,u) :- 



— 5 — . From assumption (2.5), there is a constant c > such that for 



all « > £y 2 , we have < c. Note that = f^{y 2 ,u) for u > £y 2 . So it follows 

from the assumption that we have vP~^-e~ u 



h(y,u) 2<y ? 2 u ) 



— < cw ae u g(u). By the Dominated 



Convergence Theorem, 



lim 

y^O 



u H 2 e 



Or' 



%,u) ^( 2 |^) 



du = c I u 13 a e "(in = c r(/3 - ^). 



So the claim in Proposition 2.4 follows. 



□ 
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Proposition 



2.5 



Let p = (3 — 1. For any < a < p we have that (3 — a > 1. 



S /9-"£(i) 



IS 



not integrable near zero because lim s _ 



lim 



oo by 



(1987) Proposition 1.3.6. 



Bingham et al. 



For any a > p we have that /3 — a < 1. There is a value 7 such that (3 



So lim s ^o aj a-«^i) / 
1.3.6. Therefore 



1 



lim 



) 



0, again from 



Bingham et al. 



a < 7 < 1. 
(1987) Proposition 



is integrable near for any a > p. 

1. The assertion for the second part follows from the 



sf>-<*t(\) 

Together we have the BG index is j3 
asymptotic implied by Proposition |2.4| □ 
Theorem 2.5\ The proof is entirely similar to the proof of the necessity part of Theo- 
2i4| First, JS Theorem III. 3. 24 implies ([2~7|), and (|2~9|). 



rem 



To prove the (2.10) and 

vM^s-M,2/) J (iyds by 



the form of the density process, replace J f y _^ Q yy'ir' (X s _{u)),y) 

Jo fy=/=o(V^( s ^ w > v) ~ l) 2 ' 7r (^s-( 6t; )) y)dyds and the rest remains the same. 

Theorem 2.7. (1) First we notice two facts. (1) On any compact interval I G 



□ 

there 

exists a constant C depending on /, such that for n > 1 (c.f. Nikiforov and Uvarov (1988) p. 54 
Eq. (28a)) 

\<p n (x)\ < Cn- 1/4 , x € I. (B.4) 
(2) \f n \ < 11/11 for all n by the Cauchy-Schwartz inequality. 



e Kn *n 4 which 



For t > 0, on one hand, the RHS of pTlj ) is bounded by |/o| + C||/|| ££°=i 
is finite due to the rapid decay of e~ Knt . This expansion converges absolutely for each x and 
uniformly in x on compacts, thus it defines a continuous function. On the other hand, the 
function Vtf(x) is infinitely differentiable in x. In fact, if x is replaced with a complex variable 
z, Vtf(z) is an entire function (see Theorem 3.1 in Thangavelu (2006)). The L 2 -convergence 
implies convergence almost everywhere in this case. To be more precise, let S(x) denotes the 



RHS of (2.11), and S n (x) its n-th partial sum. Convergence of S n (x) to Vtf(x) in L implies 



that there is a subsequence Sk n (x) converging to Vtf(x) almost everywhere. But the limit of 
Sk n (x) is S(x), so S(x) = Vtf(x) almost everywhere. Furthermore, since both sides of (2.11) 



are continuous functions, they must agree at every point. Therefore, for the OU semigroup, the 
eigenfunction expansion in (2.11) is valid pointwise for each / G L 2 (R,m) and t > 0. 



(2) For the SubOU semigroup, when the eigenfunction expansion on the RHS of (2.13) 



converges absolutely for each x, the spectral representation (2.13) for Vff(x) is valid for each 
x , as the following calculation can be justified: 



Vff(x) 



V s f{x)q t {ds)= / y^e- Kns f n y n {x)q t {ds) 

[0,oo) J[0,oo) n=Q 



00 „ 00 

/ e-™ s q t (ds)f n <p n (x) = e-K^f, 

n=0 ^[0,oo) n=Q 



,<Pn{X)- 



In the above, we first use the definition of the SubOU semigroup, then represent V s f{x) by the 
eigenfunction expansion which also converges pointwise, interchange the summation and expec- 
tation justified by the absolute convergence of the expansion and the dominated convergence 
theorem, and use the Laplace transform of the convolution semigroup. 



For t > 0, either condition (i) or (ii) in Theorem 2.7 ensures the absolute convergence of the 



expansion for each x, and thus the eigenfunction expansion for the SubOU semigroup converges 
pointwise. □ 
Theorem 2.8. Convergence of the expansion in the RHS of (2.14) to the RHS of (2.1) 



follows from the well-known Mehler formula for Hermite polynomials (e.g., Thangavelu (2006) 
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Proposition 2.3). Part (2) follows from the estimate (B.4) for the eigenfunctions. 



□ 



2*2 



Theorem 2.9. We first notice that for all n and all real x, \ip n (x)\ < 1.0864e ^ 2 

Therefore we have 

2 



This bound is given in Boyd (1984) and is shown there to be tight 

-0On)t| f 
\Ji 



\Zn=M^ Kn)t fnMx)\ 

that J2n=o e 



K(x-ey 
< 1.0864e^? 2_ 



v-^oo 



■M ' 



Using \f n \ < ll/H and assuming 



Lemma 



<f>(nn)t < oq i s satisfied for all t > 0, we obtain the estimate in Theorem 2.9. 

First, note that the function e x G !?{. 



3.1 



and Theorem 



3.1 



□ 

so the 



spectral representation theorem applies and 



1 



- H n 

V2"n! 

1 



0) 



-(y- 



H n (y)dy = e 6+ ^ 



dx 
1 



a 



'2k 



where we used the identity fZo e ~ (y ~ Z) H n(y)dy = y/n(2z) n (Prudnikov et al. (|1986[) p.488 
No. 17 of 2.20.3). It can be shown by using the estimate of the eigenfunctions? B.4 ) that 



the Hermite expansion of the exponential function is absolutely convergent for each x, hence 



condition (i) in Theorem 2.7 is satisfied. The results in Theorem 3.1 are obtained by applying 

□ 

3.2 Necessity. Let (B p ,C P ,H P ) be the semimartingale characteristics of the 



(2.13) to e a 



Theorem 



SubOU process with generating tuple {tip,6p,cjp : ^p,vp). Then (B p + H, C P ,H P ) is the set 
of characteristics for X under P. Since P and Q are locally equivalent, Theorem |2.5| implies 
condition (2) and (3), and that there exists some deterministic function f3 such 



Bf(u) + H(t) = B t (u) + ja 2 [ (p s (oj) + /3(s)) ds 

Jo 



+ 



J[0,t]xR 



where P s (u) 



(■ypKp0p-^K0)~(~jpKp--jK)X s ^(u) 



{ 7 ^0}- 



Thus if 7 > 0, then H is an absolutely 



0. If 7 = 0, then H(t) = for all t. 



continuous function of time, and H(0) 

Sufficiency. If 7 = 0, then the conclusion is directly implied by Theorem 2.4 
we can first find a measure 



then using Theorem 



2.4 



If 7 > 0, 
and under P, X 



locally equivalent to 

is a SubOU process with generating tuple (ftp, Op, crp,Jp,fp). Let X c be the continuous local 
martingale part of X under P. Since H is absolutely continuous, we can define X(t) := ~^r^w^- 

Then define a measure P by = S{\ ■ X c ). This is a Radon-Nikodym density process because 
the Novikov condition is satisfied, so the stochastic exponential is a true martingale. Now P 
and P are locally equivalent. Under P, the first component of the semimartingale characteristics 
becomes B P + f Q X s jpapds = B P + H(t). Thus, X is a SubOU process with the generating 
tuple (kp, 6p, ap, 7p, vp) plus a deterministic function H(t). □ 



Theorem 3.4 Since the put payoff is bounded and the measure is Gaussian, it belongs to 



L (R, m). The expansion coefficients are computed as follows: 



(K - F(x, t, t*)) + (p n (x)m(x)dx 



/oo 
(K - F(x, t, t*))l{ x<x *yip n (x)m(x)dx. 



Kl 



{x<x*} l Pn{x)m{x)dx = —= 

V 7T 



K 



2 n n\ 



(x*-6 



H n (x)e x dx 



K 



Vvr2 n n! 



b n (w*). 
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The integral in (3.9) is given in Prudnikov et al. (1986). For the second integral, 



F(x, t, t*)l {x<x * }Vn (x)m(x)dx = / Fe~ G ^ £ e~^ T f m <p m {x)v n {x) m {x)dx 



m=0 



m=0 
1 



l {x)fn(x)m(x)dx = Fe 9+ ^ G ^ Ve"' —<(„.,. 

\/iT2 n n\ — ' 

m=0 



(sm)T ^_ 

ml 



The interchange of integration and summation is justified by the Dominated Convergence The- 
orem due to the estimate: 



M 



(f m (x)(p n (x)m(x)dx 



\<Pm(x)cp n (x)\m(x)dx ^ \\(Pm\\ ■ \\ip n \\ = 1, 



and E"=o e-^ (Km)T /m < oo. With some further simplifications we obtain (|3.7|). The integral 
in ( |3.10| ) is calculated as follows. Consider the integral Jn, m ( x ) := f-oo Hn{z)H m (z)e~ z dz. By 

the identity H n {z)H m {z) = £™" (n,m) (™) (l)2 k klH n+m _ 2k (z) ( |Prudnikov et H| Jl986] p.640 
No. 11 of 4.5.1), we have 

min(n,m) , . , . min(n,m) 



k=0 

Theorem \4-l 



fc=0 V / V / -/-oo fc=Q \ / \ / 



oo 

E[/(Y t )|y s , Z a ] = E[E[/(Xr t )|T t - T s , Y s , Z s ]\y s , Z s ] = e[^ e~^ Kn ^- T ^ f nVn {Y t 



n=0 



n=0 



n=0 



where condition (1) or (2) in Theorem 4.1 justify the interchange of summation and expectation. 

□ 

Proposition 4-1 Define Z% = Zg , where S is the inverse of T. Then Z 4 C = Zx t ■ Since the 
time change S is continuous , Z c is a dapte d to S (see Jacod (1979) Definition X.13 for adaption 



to a time change), and by Jacod (1979) Theorem X.16, Z is a continuous local mar tingale 



Jacod 



w.r.t. (Ft)t>o- Now [Y c , Z c ]t = [Xj,, Z?\t = [X c , Z\x t , where the second equality is from 
(1979) Theorem X.17. Since X and Z are independent, X c and Z are independent. Because 
the cross- variation of two independent continuous local martingale is 0, we have [X c , Z]t = 
for all t, hence [X c , Z]x t = 0, and the claim is proved. □ 
Theorem \4-4\ Conditioning on the terminal state Z% of the CIR process, we have: 



roo 

E[(K-F(Y t ,Z t ,t,t*)) + ]= / E[(K-F(Y t ,z u t,t*)) + \Z t = z t ] Pcm (t,z , 

Jo 



z t )dz t 



E 



Y,e-^ Kn)Tt p n (t,f,w*,F)Myo) 



n=0 

oo 



Z t = zt 



Z t = z t 



{ 12 E [e~ H ™ 
ln=0 

/•oo I 00 

J o U=o 



Pcm(t,z ,zt)dzt 
p n (t, t*,w*,F)ip n (y ) | pcm(t, z , z t )dz t 

nn)\zo, z t )p n {t, t*,w*,F)(p n (y ) > pcm(t, z , z t )dz t . 



33 



The interchange of expectation and summation is justified by the assumption. □ 
Proposition 5.1. Using the recursion for Hermite polynomials, for n > 1, m > 0, 



/x rx 
H n+1 (y)H m+1 (y)e~ y dy = / [2yH n (y) - 2nH n ^ 1 (y)]H m+1 (y)er y dy 
-oo J —oo 

= - / H n (y)H m+1 (y)de y - 2na n _i, m+ i(x) 

J — oo 

/X 
[H' n (y)H m+1 (y) + H n (y)H' m+1 (y)]e~ y2 dy - 2na n ^ m+1 (x) 
-oo 

/X 
2 
[2nH n _ 1 (y)H m+1 (y) + 2(m + l)H n (y)H m (y)]e~ y dy - 2raa n _i )m4 
-oo 

= 2(m + l)a n>m (x) - H n (x)H rn+ i(x)e~ x2 . 

It is easy to verify that this recursion is also true for n = 0. Therefore we have 

a n+ i, m+ i(x) = 2(m + l)a n>m (x) - H n (x)H m+1 (x)e~ x2 , n>0,m>0. (B.5) 
In particular, a n ^ n (x) = 2na n -\^ n -i{x) — H n ^\{x)H n {x)e~ x , n > 1. Noting the symmetry 



fln/m^) = a m>n (x), we also obtain the following by exchanging the role of n and m in (B.5): 

a m+ i jn+ i(x) = 2(n + l)a n , m (x) - H m {x)H n+1 (x)e~ x \ m>0,n>0. (B.6) 
If m ^ n, subtracting ( |B.6 ) from (B.5), we obtain: 



a„, m (s) = e x (i? n (x)i? m+ i(x) - fl m (x)fl n+ i(a;))/(2(m - n)) (n / m, n > 0, m > 0). □ 
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